Model for the WaterâAmorphous Silica Interface: The Undissociated

J. Phys. Chem. B 2007, 111, 11181-11193

11181

Model for the Water-Amorphous Silica Interface: The Undissociated Surface Ali A. Hassanali† and Sherwin J. Singer*,†,‡ Biophysics Program and Department of Chemistry, Ohio State UniVersity, Columbus, Ohio 43210 ReceiVed: May 15, 2006; In Final Form: June 29, 2007

The physical and chemical properties of the amorphous silica-water interface are of crucial importance for a fundamental understanding of electrochemical and electrokinetic phenomena, and for various applications including chromatography, sensors, metal ion extraction, and the construction of micro- and nanoscale devices. A model for the undissociated amorphous silica-water interface reported here is a step toward a practical microscopic model of this important system. We have extended the popular BKS and SPC/E models for bulk silica and water to describe the hydrated, hydroxylated amorphous silica surface. The parameters of our model were determined using ab initio quantum chemical studies on small fragments. Our model will be useful in empirical potential studies, and as a starting point for ab initio molecular dynamics calculations. At this stage, we present a model for the undissociated surface. Our calculated value for the heat of immersion, 0.3 J‚m-2, falls within the range of reported experimental values of 0.2-0.8 J‚m-2. We also study the perturbation of water properties near the silica-water interface. The disordered surface is characterized by regions that are hydrophilic and hydrophobic, depending on the statistical variations in silanol group density.

The interface between amorphous silica and water is one of the most ubiquitous and important in chemical, biochemical, and environmental systems.1,2 Besides various forms of what is commonly known as “glass”, amorphous silica is also present as an oxidized external layer in silicon-based devices.3,4 The strong affinity of biomolecules for hydroxylated silica surfaces, and the variety of chemical modifications possible, make silica a key material for chromatographic applications.5 The strong interactions of aqueous solvent and adsorbates with the surface of silica make this system interesting and challenging. Specifically, the surface hydroxyl groups on silica can form strong hydrogen-bond complexes with water, and with biological molecules, for example, via phosphate groups.6 The interaction with certain silicas can have pathological effects on the cellular system of the human body.1,7 Silica and silicon with an oxidized surface layer of amorphous silica3 are materials commonly used for fabrication of micro- and nanofluidic devices. Electrochemical and electrokinetic function of these devices depends on the microscopic properties of the silica-water interface. Modeling crystalline and amorphous silica has been approached with empirical potentials,8-17 and also by ab initio simulation methods.18-26 Since the structure of silica is determined by making and breaking covalent bonds, ab initio methods are preferred but not always feasible. In particular, generating starting configurations for simulations of amorphous silica involves a lengthy annealing process27 which is not feasible using purely ab initio methods. As a result, all ab initio simulations of amorphous silica to date are based on starting configurations generated by an annealing process involving an empirical potential. The most popular choice of empirical potential for silica is the van Beest-Kramer-van Santen14 (BKS) model. This model correctly reproduces the structural features of many crystalline phases of silica and amorphous silica. While it might have shortcomings, especially for sur† ‡

Biophysics Program. Department of Chemistry.

TABLE 1: Pairwise Potential Buckingham Parameters atom pair

A (eV)

F (Å)

C6 (eV‚Å6)

Si-OS OS-OS OS-H OS-HW O-OW H-OW Si-OW

18003.75 1388.773 10907.00 70.795 2840.780 70.795 1049.880

0.2052048 0.3623188 0.1253934 0.3062000 0.318112 0.2662000 0.4000000

133.5381 175.0000 4.088270 0.000000 230.9240 0.000000 0.000000

TABLE 2: Three-Body Potential Parametersa atom triplet

k (eV)

θ0 (deg)

a

F (Å)

D (eV)

Si-OS-H H-OS-H OS-H-OS Si-O-Si

0.164816 13.00000 13.00000 30.00000

64.2835 0.00 0.00 0.00

-0.423721 0.00 0.00 0.00

1.85 1.4 1.4 1.6

1.24-1.5 0 0 0

a

The Si-O-H parameters control the truncated Vessal angle potential of eq 1, and the H-O-H, O-H-O, and Si-O-Si parameters are for the blocking potentials of eqs 2-4.

faces,17,28 the BKS model has been employed and evaluated in a large number of studies. Our objective is the development of a potential model that is tractable for large-scale simulations, that can take advantage of well-tested models for bulk silica and water, and that provides a pathway for integration with common force fields for biological simulation. To accomplish these goals, we have constructed our model as an extension of the BKS model14 for amorphous silica and the SPC/E model29 for liquid water. The BKS is a two-body potential and the SPC/E model is a rigid, nonpolarizable water model. Both are computationally efficient and have been used extensively. Feuston and Garofalini have previously constructed a remarkable model for the water-amorphous silica interface.13,30-33 In their model, silicon, oxygen, and hydrogen atoms are not preassigned to molecular units such as water molecules. In their model, hydrogen and oxygen atoms may evolve to be part of

10.1021/jp062971s CCC: $37.00 © 2007 American Chemical Society Published on Web 09/06/2007

11182 J. Phys. Chem. B, Vol. 111, No. 38, 2007

Hassanali and Singer

Figure 1. Ab initio calculations were performed for the single-silanol (left) and geminal (right) clusters shown here to determine parameters for our empirical potential. A subset of the capping hydrogens is indicated.

Figure 2. Fitted vs ab initio energies for (a) OH stretches of a geminal fragment and (b) SiOH angle. In the left-hand plot, points that agree perfectly would lie on the line with unit slope.

either a silanol group or water molecule, depending on the course of a simulation. For example, Feuston and Garofalini formed a hydroxylated silica surface in their simulations by placing water near a freshly cleaved surface, allowing surface chemistry to happen at 1000 K, and then annealing to ambient conditions. If the model contains all of the right “chemistry” that governs silica-water reactions, this is an attractive route to modeling the silica-water interface. However, it is not clear that the right chemistry is built into the Feuston-Garofalini model. For example, when used to model pure water, 1.4% of the water molecules engage in symmetric hydrogen bonds in which a hydrogen atom is roughly equidistant to two oxygen atoms.32 Rustad and Hay34,35 developed an improved model in this vein which was capable of yielding reasonable acid dissociation energies for orthosilicic acid in the gas and solution phase. The adsorption sites of single water molecules on the dehydroxylated amorphous silica surface have been considered by Bakaev and Steele,36 and those on the hydroxylated surface, by Leed and Pantano.37 Since we are ultimately interested in properties like diffusion and nonequilibrium flow near silica surfaces, we require

that the bulk water region be described by a model for which properties like the density, diffusion constant, and viscosity have been systematically tuned and benchmarked. Another goal of ours is to incorporate biomolecules near the silica-water interface. Therefore, instead of the FeustonGarofalini or Rustad-Hay models, more standard water models are desired because force fields for biological simulations have been developed. Unfortunately, not enough atomic-scale experimental information is known about the silica-water interface to provide a rigorous benchmark of this aspect for any model. The goal of this work is to devise a model for the hydrated, hydroxylated silica surface as an appropriate extension of the BKS model for amorphous silica14 and the SPC/E model for water.29 To that end, we have attempted to cast the interaction model as much as possible in a form compatible with the BKS and SPC/E models. However, surface behavior is distinct from bulk behavior, and we have found it expedient to include some three-body interaction terms that are not found in the BKS potential. Since it is not possible to predict ahead of the surface annealing process which silicon and oxygen atoms will be part

Model for the Water-Amorphous Silica Interface

J. Phys. Chem. B, Vol. 111, No. 38, 2007 11183 1. Development of a Model for Hydrated Amorphous Silica

Figure 3. Configuration in which a hydrogen atom is simultaneously within covalent bonding distance of two oxygens. Such divalent hydrogen configurations, as well as spurious -SiOH2 groups, are eliminated by the HOH and OHO blocking potentials, eqs 2 and 3.

of surface silanols and siloxane bridges, and because we aimed for a simple potential model with fixed partial charges, we constrained the charges of all silicon and oxygen atoms to be those found in the BKS model. The parameters of our model were determined using ab initio quantum chemical studies on small fragments. In this respect, the procedure used to develop our model is opposite of the one employed by Cruz-Chu et al.,38 who optimized interaction parameters by matching a bulk property, the contact angle of a water drop. This procedure could only be carried out when the surface was sufficiently hydrophobic to have a nonzero contact angle which, in practice, did not allow optimization when the surface contained silanol groups. Our model will be useful in empirical potential studies, and as a starting point for ab initio calculations. The model presented here should be regarded as an initial version to be calibrated and refined by comparison with ab initio simulations and experiment. The dissociation of silanol groups is an important chemical feature of the hydroxylated silica surface, leading to its characteristic negative charge. This requires a significant extension of our potential, which is currently under development. At this stage, we present a model for the undissociated surface. In section 1, our procedure for developing our model is described in detail. The functional form of the potential and ultimate parameter values are given. This functional form should be regarded as a minimal form required to obtain reasonable physical behavior. In section 2, results of simulation of the silica-water interface are presented. The protocol for generation of the slab is described. The heat of immersion is an important benchmark. Our model is shown to fall within the range of available experimental values. We also report structural features of the silica-water interface that result from our interaction model and protocol for hydroxylating and annealing the silica surface. Finally, we discuss the implications of our work and prospects for future development.

1.1. Formulation of the Potential. So we could rely on a well-tested model of the bulk, we designed our model of the hydrated silica surface as an extension of the BKS14 and SPC/ E29 potentials. The BKS potential is a sum of pairwise Coulomb and Buckingham (exponential repulsion + r-6 attraction) interactions. In the BKS model, silicon and oxygen atoms are assigned charges of +2.4 and -1.2, respectively. To maintain a model with fixed charges, we kept the same charges for surface atoms, even though in reality the BKS charges are artificially large and, in any case, the effective partial charges at the surface are expected to be different from the bulk. In order to avoid overbinding of water to the silica surface, we found it necessary to distinguish between silanol type and siloxane type oxygens, for the short range potentials, while using the same charge of -1.2. As shown below, an acceptable fit to ab initio data could be achieved with these fixed charges. Charge neutrality of the undissociated surface requires that the charge on hydrogen atoms be +0.6. It is remarkable that the BKS potential describes many crystalline forms as well as amorphous SiO2 using no more than pair interactions.14 We were not successful devising a potential for the surface based on pair interactions, and found it necessary to include three-body terms. The most important is a SiOH angle bending term, which we found to be described well by the truncated Vessal form:39,40 uTV(θ, rSiO, rOH) )

a k θa(θ - θ0)2(θ + θ0 - 2π)2 - πa-1(θ - θ0)2(π - θ0)3 - DTV × 2 exp[-F-8(rSiO8 + rOH8)] (1)

[

]

where θ is the Si-O-H angle, rSiO and rOH are the two bond lengths, and there are three adjustable parameters k, a, and θ0. The usual truncated Vessal potential is zero at θ0 and positive elsewhere. The exponential factor, which cuts the interaction at any θ * θ0 from a positive value to zero when either rSiO or rOH becomes large, thereby provides a large energetic incentive for the SiO bond and especially the OH bond to break. While the OH bonds are in principle dissociable, the current form of our model does not encompass dissociative behavior, least of all through a cutoff in an angle potential. Therefore, a negative energy contribution of -DTV is added to balance the reward for dissociation that arises as an artifact because the otherwise positive truncated Vessal potential is cut off at large separations. While a combination of Coulomb, Buckingham, and SiO-H bending potentials gave an acceptable fit to ab initio data (see below), we found that thermal simulations of a model with just these interactions led to unphysical associations, hydrogen atoms that were simultaneously strongly bound to two oxygen atoms, or silanol oxygen atoms that accepted a second hydrogen, i.e., -SiOH2 groups. We blocked out unphysical O-H-O and H-O-H association using two more threebody terms given in eqs 2 and 3. Also included is a Si-O-Si blocking potential that prevents silanol oxygens from forming three coordinated species on the surface shown in eq 4.

uOSHOS(rOSH, rO′SH) ) k exp[-F-4(rOSH4 + rO′SH4)]

(2)

uHOSH(rOSH, rOSH′) ) k exp[-F-4(rOSH4 + rOSH′4)]

(3)


uSiOSi(rSiO, rOSi′) ) k exp[-F-4(rOSi4 + rOSi′4)]

(4)

In the formulation of our potential, the symbol O with a subscript “X” refers to a siloxane-bulk type oxygen, the symbol O without a subscript refers exclusively to silanol oxygens, and O with a subscript “S” refers to a general silica oxygen, either silanol or siloxane. Hence, eqs 2 and 3 apply to all silica oxygens, while eq 4 applies only to silanol oxygens. The total hydroxylated silica interaction potential is of the following form.

Usilica ) UBKS + Usurface

UBKS )

NSi i-1 q 2 Si

∑ ∑ i)2 j)1 r

NOS i-1

+

ij

∑ ∑ i)2 j)1

[

qOS2 rij

∑ ∑ i)1 j)1 Usurface )

∑ ∑ i)1 j)1

[

q O qH rij

rij

]

+ uB(SiOS)(rij) (6)

]

NOS i-1 NH

uTV(θijk, rij, rjk) +

∑ ∑∑uO HO (rik, rjk) + i)2 j)1 k)1 S

NH i-1 NOS

∑ ∑∑ i)2 j)1 k)1

[

qSiqOS

+ uB(OH)(rij) +

NSi NO NH

∑ ∑∑ i)2 j)1 k)1

]

+ uB(OSOS)(rij) +

NSi NOS

NO N H

(5)

S

NSi i-1 NO

uHOSH(rik, rjk) +

∑ ∑∑uSiOSi(rik, rjk) i)2 j)1 k)1

(7)

In eqs 6 and 7, uB(r) is the Buckingham potential, a combination of exponential repulsion and r-6 attraction. Since the r-6 term diverges as r f 0, an unmodified Buckingham potential has an unphysical maximum at small r, and then plunges to -∞ at even smaller r. Therefore, it is standard36,41 to modify the Buckingham potential at short range to eliminate the unphysical maximum and make it smoothly repulsive at small r. For the interactions in Table 1 where C6 is not zero, we implemented a small r patch of the form

uB(r) )

{

Ae-r/F -

C6 , r g r0 r6

C12 , r12

r < r0

B+

}

(8)

where the switching radius, r0, is 1.1 times the location of the maximum of the unmodified Buckingham and the parameters B and C12 are chosen to the value and derivative of the potential continuous at r0. The water-silica interactions in our model are also sums of Buckingham interactions. In practice, we found that we could dispense with the London dispersion term between the silanol hydrogens and water molecules, and eliminate the Buckingham potential altogether for interactions of water hydrogen atoms with either silicon or silanol hydrogen atoms because a large Coulomb repulsion already precludes close contacts.

[

NSi NOW

Usilica-water )

[

NOS NHW

∑ ∑ i)1 j)1

∑ ∑ i)1 j)1

qOSqHW rij


qSiqOW rij

] ∑∑[

+ uB(OSHW)(rij) +

] ∑∑[

uB(HOW)(rij) +

NO NOW

qO qO W

i)1 j)1

rij

NSi NHWqSiqH

∑ ∑ i)1 j)1

]

+ uB(SiOW)(rij) +

rij

W

NH NHWq

+

∑ ∑ i)1 j)1

NH NOW

qHqOW

i)1 j)1

rij

+

]

+ uB(OOW)(rij) +

H qH W

rij

+

NOXNOWq q OX OW

∑ ∑ i)1 j)1

rij

(9)

In the above equations, the subscript “W” indicates atoms that are part of water molecules. We took the interaction between water molecules as the SPC/E model29 and used the SPC/E charges in the interaction potential. In principle, other waterwater potentials could be used in place of SPC/E. However, using different charges on the water molecules would presumably degrade the quality of fit to ab initio water-silica interactions. Alternatively, one could introduce an inconsistency between water-water and water-silica Coulomb interactions. The procedure used to fit the potentials and the quality of fits are described below. Here, for the convenience of those wishing to use the potentials, the final parameters are collected in Tables 1 and 2. As before, O with a subscript “X” refers to oxygen types that are not connected to hydrogens, O without a subscript refers to silanol oxygens, and O with a subscript “S” includes all silanol and siloxane oxygens. Besides the Coulomb potential, there are no other interactions between pairs not listed in Table 1, or triples not listed in Table 2. 1.2. Adjustment to Match ab initio Data. The free parameters in eqs 5-9 were adjusted to best match ab initio quantum mechanical calculations for the single-silanol and geminal fragments shown in Figure 1, with or without a neighboring water molecule. These clusters were excised from a silica slab taken from molecular dynamics calculations. The dangling oxygen atoms that would be connected to the rest of the slab were capped with hydrogen atoms. This section is devoted to describing the quality of that match. The BornOppenheimer energy surface was obtained using MP2 perturbation theory42-46 to account for electron correlation. The 6-311G** basis set was adequate to describe the energetics of silica surface fragments, and energy differences were negligibly changed when the basis was improved to the 6-311++G** level. However, potentials involving interactions between water molecules and silica fragments required the larger 6-311++G** basis set. All electronic structure calculations were performed using Gaussian 03.47 One of the challenges of developing robust empirical potentials is ensuring an adequate sampling of the configurational space that will eventually be explored in statistical simulations. Our initial attempts to use configurations from molecular dynamics (MD) simulations in a force matching scheme48 were not successful for two reasons. First, our initial parameter guesses generated configurations that were too far from the ones that eventually proved to be most probable, and the desired configurations were not sampled adequately in MD simulations. We then sampled a series of configurations on coordinate grids to generate better initial parameter guesses. Even with those better initial parameters, iterating the parameters to match the energies of configurations sampled from MD simulations did not lead to physically reasonable potential functions. We eventually concluded that coordinate grids were

Model for the Water-Amorphous Silica Interface most effective for fitting the type of potential functions proposed in section 1.1. Tabacchi and co-workers49 report similar difficulties in their efforts to fit potentials for NaCl by force matching. They also conclude that an extensive force database is required to develop physically meaningful parameters for their empirical potentials. Like Tabacchi and co-workers, we also attempted to expand the configurational search space by raising the temperature of the system, but this did not aid in getting physical potentials. While Izvekov and co-workers50 have been successful in using the force matching technique, their work involved potentials that depend linearly on parameters which is not true for most of our model potentials. All of the configurations that were used to fit our potentials were generated from grids that sampled a broad region of configuration space. 1.3. Silanol Groups. The short range potential functions for silanol -SiOH groups consist of SiO and OH stretching potentials and a SiOH angle bending term. The SiO stretch is constrained to be that of the BKS model, facilitating use of our potential as an extension of BKS. The configurations used to generate fitting points for the silanol potentials were based on the geminal-silanol fragment shown in Figure 1. These clusters were obtained from a larger slab. All of the clusters were capped with hydrogen atoms where they would join the rest of the slab. The OH stretch and SiOH angle parameters were fit using a geminal cluster, as shown in Figure 1. The OH distance for each -SiOH group on the geminal was varied between 0.7 and 1.3 Å, forming a two-dimensional grid. The angle grid was formed by fixing the OH distance to approximately 0.95 Å which is the equilibrium bond length for the OH stretch, and the angle was varied between 1.1 and π radians. Typical results comparing ab initio and fitted energies are shown in Figure 2. In Figure 2a, the stretching potential is tested for a geminal fragment. The complicated coupling between the two OH groups is captured in an approximate way by the various Coulomb and short range terms. Figure 2b shows how the truncated Vessal form (eq 1) captures the angle bending behavior. Although the Coulomb, bond stretching, and angle bending potentials described the most stable regions of the silanol fragment (Figure 2), this empirical potential gave rise to other, spurious minima that did not correspond to actual physically stable configurations. MD simulations performed with only the terms in the potential mentioned above gave rise to configurations in which surface hydrogens divalently bonded to two oxygens, particularly in geminal silanols because of proximity effects. This is shown in Figure 3. However, the divalent bonding also occurred between two silanol groups that were not members of the same geminal pair. This effect is not surprising given that the OH interaction is strong. We also discovered that certain oxygens could acquire a second hydrogen, resulting in a -SiOH2 group. In order to eliminate these unphysical minima on the potential surface, we introduced “blocking” potentials: HOH and OHO three-body terms which raise the energy of divalent hydrogen or -SiOH2 configurations. These terms (eqs 2 and 3) have no angle dependence. They become highly repulsive when a hydrogen atom is simultaneously within covalent bonding distance of two oxygens, or when an oxygen is within covalent bonding distance of two hydrogens. The exponential range parameter, F, is set to a value that decreases these blocking potentials to nearly zero when at least one of the oxygen-hydrogen distances is typical of noncovalent hydrogen bonds. The blocking potentials were not part of a fit to ab initio data. Instead, the blocking potential parameters were adjusted by trial and error until unphysical

J. Phys. Chem. B, Vol. 111, No. 38, 2007 11185 configurations were no longer observed and physical configurations were weakly affected. There is considerable room for adjustment of the blocking potential parameters without significantly affecting the properties of the model. Moreover, in simulations, the blocking potentials are most important during equilibration or annealing at high temperature. The additive constant, DTV, in eq 1, and the three-body blocking potentials in eqs 2 and 3 are less important for simulations near 300 K and might be reduced after equilibration and annealing depending on the nature of the surface. In our simulations, we also found the formation of a hydrogenated siloxane species where the silanol group is shared between two silicons. Ma and coworkers have reported the formation of this species in their ab initio simulations but also found that on addition of water the hydrogenated species went away, resulting in the formation of a silanol group.51 It is still unclear whether this species will exist on a realistic surface, but given the scarcity of evidence for this species in the literature, we discourage formation of this species in our model with the Si-O-Si three-body potential. Only 5% of silanol oxygens are in close proximity to a second silicon atom. 1.4. Water-Hydroxylated Silica Interactions. The watersilica potentials were fit using various possible paths of water approach to the isolated silanol group shown in Figure 1. By choosing likely paths of approach, we fix important regions of the potential surface. Initial testing reported below suggests that the potential provides a reasonable interpolation between the configurations where it was matched to ab initio data. In choosing paths of approach, we followed Saengswang et al.52 and used (1) silanol donor-water acceptor, (2) water donorsilanol acceptor, and (3) simultaneous donor-acceptor paths. These three paths are shown in Figure 4. The water-silanol separation was varied from 0.5 to approximately 6 Å to get a comprehensive sweep of the potential energy surface. The water model geometry was the SPC/E water model29 which was also used in our MD simulations. Further simulations and the resulting thermodynamics of wetting indicated that applying interactions determined from water-silanol group interactions (as in Figure 4) to all oxygens on the surface leads to overbinding of water to siloxane oxygens on the surface. We thus distinguished between silanol oxygens and siloxane oxygens. Shown in Figure 5 are three paths of approach to a siloxane cluster which we used in ab initio calculations to calibrate water-siloxane interactions. Notice that in these clusters the silicons are capped with hydrogens. In the entire set of water-hydroxylated silica pair interactions, only the water oxygen-silanol oxygen (O-OW in Table 1) interaction carries a London dispersion force. As in the construction of most water potentials, the dispersion forces centered on the hydrogen atoms are neglected. Furthermore, in typical configurations, water oxygens are not close to silicon atoms. Therefore, we also neglected a dispersion interaction between silicon atoms and water. It was still necessary to retain the exponential repulsion between water oxygens and all atoms of hydroxylated silica and between silica oxygen and water hydrogens. Without the exponential repulsion present in H-OW, Si-OW, and O-HW interactions, the attractive Coulomb infinity is discovered within a few picoseconds of MD simulation, leading to unphysical dissolution of hydrogen, silicon, and oxygen atoms from the slab. Even with the exponential repulsion, the Coulomb singularity still exists. However, the exponential repulsion causes the system to avoid the region of the attractive Coulomb singularities. We did not find it necessary to block out this singularity with a short range patch as in eq 8,



Figure 4. Paths of approach for water near silanols.

Figure 5. Paths of approach for water near siloxanes.

although this certainly could be done if needed, say, for simulations at high temperatures. Although not designed for a realistic description of aqueous dissociation chemistry, the potential form given in section 1.1 is in principle dissociable. Part of the freedom to dissociate Si-O bonds is accessed during annealing of the initially hydroxylated surface and leads to physically reasonable configurations. However, our potential is not designed to model the energetics of dissociation. The Buckingham + Coulomb form of our potential, together with the many constraints on partial charges (BKS charges for silica Si and O, SPC/E charges for HW and OW) we imposed to mesh our interface model with established bulk potentials, does not have enough flexibility to simultaneously give a good fit to all three paths in Figure 4. Path 1 leads to the waterhydroxylated silica complex with the highest binding energy, and is therefore a probable approach that a water molecule would make to the silica surface. Therefore, we weighted it more heavily than the other paths. Figure 6 shows all of the ab initio energies for all three paths of approach, and the predictions of our empirical potential for each of those paths. As can be seen from the fits, paths 1 and 3 give decent agreement with ab initio. Our model is not sufficiently flexible to simultaneously yield high quality fits to all four paths of approach. We tested the transferability of our potential model, parametrized for approach of water to a single silanol, by comparing the empirical potential and ab initio results for the approach of a water molecule to a geminal fragment. Ferrari et al. have done a more extensive study on the energetics of water near geminal type silanols.53 Our path of approach shown in Figure 7 was adopted from their ES-2-W scheme. As can be seen in Figure 6, the interaction of water with the geminal fragment is similar to that in the single silanol, and hence, our potential gives a

Figure 6. Fitted (open symbols) vs ab initio (filled symbols) energies for the paths of approach shown in Figure 4 and for approach of a water to a geminal fragment along the path shown in Figure 7. The four paths are vertically offset from each other on the graph for clarity.

qualitatively correct representation of this interaction. In principle, there are many situations that should be checked for transferability. Testing of all of these situations is beyond the scope of this work, although to some degree the verification of macroscopic properties like the heat of immersion (section 2.2) suggests that none of the interactions are seriously in error. The limited flexibility of our empirical potential and the constraints of matching the BKS and SPC/E model charges precludes a physically meaningful fit to water interactions with capped silicate clusters in all possible situations. An example


J. Phys. Chem. B, Vol. 111, No. 38, 2007 11187

Figure 8. Fitted empirical (open symbols) and ab initio energies (filled symbols) for water approaches to siloxane fragments. Our model was designed to match the more important interactions with silanol groups (see Figure 6) and, as shown here, is not sufficiently flexible to also match ab initio data for approach of waters to this siloxane fragment. Figure 7. Path of approach of water to a geminal silanol fragment.

TABLE 3: Structural Properties of Gas Phase Orthosilicic Acid (Angles in Degrees and Distances in Å) parameter

this work

Sauer54

r(SiO) r(OH) ∠SiOH ∠OSiO ∠HOSiO

1.615 0.95 114.69 106.76 -10.88

1.63 0.95 117.15 115.79 -33.34

is the poor job reproducing the interaction potential between a water molecule and the siloxane-containing fragment shown in Figure 5. Approach of a water molecule to the siloxane group in Figure 5 is dominated by the large positive BKS model charges on the two silicon atoms. The combination of the large charges from the BKS model and SPC/E water model make it impossible to avoid overly attractive interactions when the water oxygen is turned toward the siloxane group (path 1) and overly repulsive interactions when the water hydrogen faces the siloxane (path 2). Further testing of our model in condensed phase simulations indicated that this effect is mitigated when water approaches a larger fragment with a silicon-oxygen ratio closer to 0.5 than 2. Silanol interaction parameters are much more important than siloxanes in determining the properties of the aqueous interface, as shown below in section 2, and therefore, we employed the limited flexibility of the potential model to more accurately describe water-silanol interactions. 1.5. Orthosilicic Acid. While our potential is optimized with emphasis on silanol groups on the bulk surface and their interactions with water, for completeness, we document how well the potential describes isolated orthosilicic acid because this species has been considered on several occasions in past studies. (See ref 55 and references therein.) The lowest energy structure predicted by our potential is the S4 structure also found in ab initio calculations.54 Geometrical parameters for the minimum energy gas phase structure of orthosilicic acid are given in Table 3. The bond lengths and angles compare well to previously published ab initio data,54 the largest deviations occurring for the HOSiO dihedral angle.

2. Slab Studies 2.1. Slab Generation. In order to test various structural and energetic properties predicted by our empirical potentials, we conducted a series of molecular dynamics simulations of a hydrated silica slab with a hydroxylated surface. In this section, we present the methodology used to generate the amorphous silica slab and afterward present an analysis of the structural and thermodynamic properties of the slab interacting with water. All simulations were conducted using the DLPOLY package which is particularly suited for materials simulations.56 A time step of 1 fs was used. The SHAKE algorithm was used when rigid waters were present. Electrostatics were calculated using the SPME method.57,58 A cutoff of 9.0 Å was used for the Buckingham interactions, and a cutoff of 3.0 Å, for the threebody potentials. The starting structure was bulk crystalline tridymite with no free surfaces. Its starting dimensions were 33.5 Å × 34 Å × 39 Å, consisting of 896 silicon atoms and 1792 oxygen atoms. The protocol that was used to generate amorphous silica from the starting material was adopted from cycle I-IV of Huff and co-workers.27 It began with heating to 8000 K followed by annealing at 4000 and 2000 K. We cleaved the surface by opening a gap in the z dimension and followed with annealing for a limited time at 300 K. Shown in Table 4 below are the population of nonbonded oxygens (NBOs), two-membered (2M) rings, and three-coordinated silicons on the surface as a function of annealing time following cleavage. Also shown is the silanol density that resulted from different annealing times after hydroxylation with the procedure described below. Variable annealing time subsequent to cleavage provides a mechanism for controlling the surface population of geminal, vicinal, and isolated silanols. Unfortunately, the detailed structure of the aqueous silica surface, its mechanism of hydroxylation, and the distribution of silanol groups between isolated, vicinal, and geminal silanols has not been conclusively established. Because we are not performing ab initio simulations, and realistic ab initio simulations to anneal a silica surface and react it with water are not currently feasible, we generated a hydroxylated starting configuration via two transformations of the cleaved and partially annealed surface: 2M rings were reacted to form vicinal silanols, and NBOs reacted to form geminal silanol pairs. Ceresoli et al.28 have studied 2M rings on the silica surface. Bakos and



Figure 9. Rare defects taken from the amorphous silica surface with a silanol density of 6.4 nm-2: (a) NBO; (b) three-coordinated oxygen; (c) silicon with three hydroxyl groups. Defects are indicated with arrows.

experiments by Maciel and co-workers were performed either on silica in contact with water vapor or on samples evacuated at room temperature or 200 °C.65-67 They determined that geminal silanols comprise only ∼6% of the total silanol population. In contrast, second harmonic generation (SHG)68 and evanescent wave spectroscopic experiments69 on silica in contact with liquid water indicate that ∼80% of the silanols are geminals. Also, X-ray photoelectron spectroscopy (XPS) studies of polycrystalline quartz in contact with solutions of varying pH found a nearly constant 1.8:1 ratio of surface oxygen to surface silicon atoms, indicating that most surface silanol sites were geminals.70 There is considerable justification for transforming nonbridging oxygens (NBOs) on the partially annealed cleaved surface to a geminal silanol by the reaction

Figure 10. Hydrogen-bonded silanol chain on our silica surface.

co-workers,59 Masini and Bernasconi,60 Du et al.,61 and Mischler et al.62 have performed ab initio and mixed quantum-classical simulations to study how water interacts with amorphous silica in the bulk and have found that water interacts with 2M rings on the silica surface to generate two close-by silanol groups, known as vicinal silanols. Following the mechanism elucidated in these works, we converted each 2M ring into vicinal silanols according to the reaction

For lack of more detailed information on the mechanism of hydroxylation, we use the NBOs to generate geminal silanols. Reviewing the situation in 1979, Iler1 stated that geminal silanols probably do not exist on a dried silica surface.63 Despite his skepticism about geminal silanols on dried silica, in aqueous solution, Iler raises the possibility that aqueous solution chemistry can raise the number of hydroxyl groups bound to a surface silicon atom.64 Subsequent NMR experiments indicated that geminal silanols are present on the surface of both dried silica, and silica in contact with water vapor.65-67 The NMR

This procedure was also followed by Rignanese et al.71 It is likely an oversimplification. Walsh et al. have investigated the reaction of NBOs with water, finding that the reaction can lead to either geminal or vicinal silanols, depending on the surface structure adjacent to the NBO.72 Also, the properties of nonbridging (dangling) oxygen atoms have been investigated by Murashov,73 who illustrates that they are not a uniform species. Our hydroxylation and annealing procedure, as described below, brings several properties of our surface in line with experimental measurements, but it must be regarded as preliminary pending further theoretical and experimental efforts to characterize the hydroxylated amorphous silica surface. The actual chemistry, as indicated above, is almost certainly more complex. Significant reconstruction of the surface took place during annealing of the initial hydroxylated configuration described above. The initial hydroxylated surface was quenched to a local potential energy minimum before further annealing. The quenching process removed some particularly high energy and presumably unrealistic features. For example, in rare instances, the same silicon atom has an NBO attached to it and is also part of a 2M ring. Under our hydroxylation protocol, this leads to a presumably bad starting configuration with three hydroxyl groups attached to one silicon. We have found that quenching of these high energy configurations can induce reconstruction of these over-hydroxylated silicons to more reasonable configurations. Alternatively, excess silanol groups can be transferred to three-coordinated silicons. Ab initio simulations by Ma et al.51 provide some evidence that such reconstruction is realistic. Similarly, geminal silanols can also


J. Phys. Chem. B, Vol. 111, No. 38, 2007 11189

Figure 11. Number densities for water and siloxane oxygens (solid curve and long dash curves, referred to left axis) and silanol oxygens (short dash curve, referred to right axis) along a direction perpendicular to a silica slab with a surface silanol density of 6.4 nm-2.

Figure 12. Figure showing affinity of water for hydrophilic regions on the silica surface. See discussion in text.

form on silicons that are five-coordinated, and some care may have to be taken in the initial few picoseconds of the simulations. Adjusting the magnitude of DTV to a more negative value is useful at this stage to ensure that the energy made available by relaxation of the surface does not induce dissociation of fragments from the surface. Once hydroxylated, the surface was resistant to further “chemical” change. We analyzed hydroxylated surfaces for unusual bonding configurations and found they were quite rare. For example, the surface with a silanol density of 6.4 nm-2 exhibited three NBOs (Figure 9a) that reformed after all of the NBOs following cleavage were reacted according to eq 11. Also, we found three three-coordinated oxygens (Figure 9b) and a single silicon that was triply hydroxylated (Figure 9c). The vast majority of the bonding maintained its original pattern. Due to the scarcity of resources in the literature, the actual concentration of these surface defects remains an open question.

For our heat of immersion calculations, we generated two surfaces with different silanol densities using the protocol described earlier. The surface with a higher silanol density of 6.4 nm-2 is in good agreement with those quoted in some previous studies.1,74 This surface yields ∼70% geminals and ∼30% vicinal silanol groups. Surface relaxation can lead to isolated silanols as well. Because the silicons prefer to be four-coordinated, surface relaxation due to an overcoordinated silicon can result in the transfer of a hydroxyl group to an undercoordinated silicon but only occurs to our knowledge at the beginning of the simulations when the surface is at a high potential energy. This is in good agreement with second harmonic generation (SHG)68 and evanescent wave spectroscopic experiments69 but with a much higher fraction of geminals than that found by NMR spectroscopy.65-67 Rarivomanantsoa et al. have found a significant number of NBOs at a surface governed by the BKS potential.75 They report that


Hassanali and Singer TABLE 4: Population of NBOs, 2M Rings, and Three-Coordinated Silicons on a Silica Surface as a Function of Annealing Time of a Freshly Cleaved Surface at 300 K annealing time (ps)

NBOs

2M rings

three-coordinated silicons

silanol density (nm-2)

0.10 0.25 1.65 2.00

48 28 17 16

24 32 31 26

18 13 4 4

6.60 5.56 4.45 3.89

TABLE 5: Average Potential Energy for Systems Consisting of 896 Silicon Atoms and 3220 Water Molecules Used to Calculate the Heat of Immersiona Figure 13. Radial density of water molecules surrounding silanol oxygens (solid curve) and surface siloxane oxygens (dashed curve). Siloxane oxygens were identified as surface siloxanes if they were at least 12.5 Å from the center of the slab. (See Figure 11.) Asymptotically, the density of water oxygens is approaching one-half the bulk water density ((1/2)0.0335 Å-3) because there are no water molecules in the half space below the surface.

15% of the surface oxygens are NBOs. The surface with the lower silanol density of 4.0 nm-2 also falls into the range of silanol densities quoted by other workers. This surface yields 26 2M rings on the surface and 17 NBOs. The relative populations of geminals and vicinals are now ∼40 and ∼60%, respectively, in better agreement with NMR studies, as mentioned earlier. Clearly, more theoretical and experimental work is needed to better characterize the silica surface. Besides resolving inconsistencies in experimental data, better theoretical annealing of the silica surface, both with respect to better annealing protocols and extending treatment to ab initio methods, is needed. For example, Ceresoli and co-workers28 find that performing the annealing using the BKS potential generates good bulk silica but generates a surface that has too many overand undercoordinated silicons and oxygens. However, these surface properties depend upon annealing protocol and sampling temperature.75,76 The annealing procedure for the cleaved silica surface affects the starting point for our hydroxylation procedure, and ultimately the characteristics of the hydroxylated surface. 2.2. Heat of Immersion. Since the properties of amorphous silica depend on preparation, it is not surprising that the heat of immersion varies between different experiments and different samples.1 Measurements of the heat of immersion date back to the work of Boyd and Harkins, who reported a value of 0.6 J‚m-2. Taylor and co-workers77,78 report the heat of immersion as 0.16 J‚m-2 and claim that it is independent of the sample specific surface area. However, Makrides and Hackerman79 demonstrate a specific surface area dependence of the heat of immersion and quote ranges between 0.16 and 0.88 J‚m-2 depending on the type of dehydration pretreatment to rid the system off extraneous water. In fact, Makrides and co-workers quote a value as high as 1.45 J‚m-2 for one of their surfaces. Chikazawa and co-workers80,81 report heat of immersion values for various silicas in the range 0.157-0.23 J‚m-2 and suggest that this value is quite sensitive to the degree of crystallinity of the surface as well as the density of surface silanols. Clearly, the experimental values cover a range of values, due to differences in porosity and specific surface area, different degassing temperatures, and the type and extent of defects on the surface. Despite the variability in experimental values, the heat of immersion provides an important benchmark for our model.

system

〈energy〉 (4.0 nm-2)

〈energy〉 (6.4 nm-2)

dry silica slab (eV) water (eV) hydrated slab (eV) enthalpy change (eV) heat of immersion (J‚m-2)

-52397.38(0.084) -1596.80(0.12) -54031.81(0.62) -37.63(0.82) 0.279(0.0061)

-52962.60(0.073) -1596.80(0.12) -54671.62(0.42) -112.22(0.62) 0.832(0.0046)

a Error bars for each run are estimated by the blocking method82 by partitioning each run into four blocks. The combined surface area of the two sides of the silica slab is 2.16 × 10-17 m2. The two sets of data are for surfaces containing silanol densities of 4.0 and 6.4 nm-2.

The theoretical calculation for the heat of immersion requires three separate simulations: dry silica slab, pure water, and hydrated slab. The dry slab simulation was done with NVT dynamics and a separation of 60 Å between the top of one slab and the bottom of its periodic replica in the next cell. The pure water and hydrated slab simulations were conducted using NPT dynamics. They contained 3320 waters. For the slab-water simulation, we generated a starting condition by inserting waters into the system no closer than 3.5 Å from the slab. During the subsequent NPT dynamics, the water approached the slab and penetrated it to some extent, as described in the following section. Water did not penetrate to the voids deep in the bulk that BKS model silica is known to contain. These voids were not filled by other means in the hydrated slab simulation. This is justified because such deeply buried water would not be removed during the drying process and would not contribute to the heat of immersion. The pressure-volume contribution to the enthalpy of immersion was calculated and found to be negligible compared to the energy. The data used in the calculation of the heat of immersion are given in Table 5. In the table, we show the heat of immersion for two surfaces with different silanol densities. The surface with the lower silanol density yields a heat of immersion that is in better agreement with most of the reported experimental values. However, the heat of immersion for the higher silanol density system also falls at the higher range of quoted heats of immersion. Long range corrections for the r-6 dispersion energy, which are only approximately implemented for inhomogeneous systems, were a minor correction to the total energy, and were not included in the numbers reported in Table 5. 2.3. Structural Features of the Hydroxylated Silica-Water Interface. The radial distribution functions (rdf’s) for silanol O-H and Si-O separations indicate that the equilibrium O-H and Si-O distances are, as expected, approximately 0.96 and 1.63 Å, respectively. The Si-O rdf essentially reflects the contributions of bulk Si-O pairs. However, these distances do not significantly change on the silica surface. The Si-O-H angles vary depending on the surface morphology and type of silanol. Our surface yields a fairly broad range of angles between 105 and 130°. Previous simulations and electronic structure

Model for the Water-Amorphous Silica Interface calculations yield a Si-O-H angle of about 121°. The smaller angles sampled on our surface are indicative of pockets of silanol clusters on the surface that form linear hydrogen bonds, as shown in Figure 10, which effectively reduces the Si-O-H angle. Shown in Figure 11 are the number density distributions for silanol, siloxane, and silanol oxygens along a direction perpendicular to the surface of a slab with a silanol density of 6.4 nm-2. Our z-densities for water relative to the surface are somewhat sharper than those calculated by Warne et al.83 Regions on the silica surface with high silanol densities and hence closely knit hydrogen-bonded networks (Figure 10) serve as regions of higher water density and potentially higher heats of immersion. This is consistent with results from Takei and Chikazawa80 and Fubini et al.84 who propose that hydrogenbonded clusters of silanols can lead to higher heats of immersion due to stronger water-silica interactions. This will be confirmed in our discussion of Figures 12 and 13. Further evidence for the affinity of water to regions of higher silanol density is seen in Figure 12 which illustrates the existence of hydrophilic regions and hydrophobic patches on our silica surface. It should be noted that Fubini and co-workers85 distinguish hydrophobicity and hydrophilicity of the silica surface according to the presence of different silanol groups (vicinal, isolated) on the surface. While this may have an effect on the binding energy of water to the specific silanol, we distinguish the hydrophobic and hydrophilic regions on the basis of overall silanol density. Nevertheless, we demonstrate the heterogeneity of the silica surface, as do they. Figure 12 is a snapshot looking down into the slab in the z-direction. The patches outlined in the left panel of Figure 12 contain surface siloxane bonds only with no silanols, and are apparently hydrophobic, as seen in the distribution of surface waters shown in the right-hand panel. This effect is quantified in Figure 13 which depicts the difference in the radial density of water oxygens surrounding silanol oxygens and surface siloxane oxygens. Furthermore, in Figure 12, one can also see that the presence of voids can lead to water clustering and possible channels for water trickling and diffusing into the surface. Bakos and co-workers59 have elucidated mechanisms for water diffusion through amorphous silica. 3. Conclusions Further understanding of the surface properties of amorphous silica is central to the theory of electrochemistry and electrokinetic phenomena, and important for a variety of applications. The silica surface figures prominently in the literature on the electric double layer and induced electrokinetic effects.86,87 The availability of reactive silanol groups furnishes a strategy for tailoring the silica surface by chemical modification for selective metal ion extraction.88 The surface potential of silica leads to a electroosmotic flow,87 which is a useful way to induce fluid flow in micro- and nanochannels.89 The surface potential can be controlled by pH and chemical modification. Covalent attachment or physical adsorption of biomolecules to silica is a common strategy for analytical and biomedical applications.90-97 In the case of human carbonic anhydrase, the regions of the protein that bind to silica nanoparticles have been established.98 The effect of succinylation on the adsorption of lysozyme to the amorphous silica surface of an oxidized silica wafer has been reported.99 Time-resolved fluorescence has been used to monitor the adsorption on silica particles of small peptides,100 including the Lys-Trp-Lys peptide, a target of our recent theoretical study.101

J. Phys. Chem. B, Vol. 111, No. 38, 2007 11191 The toxicity of silica, which is highly dependent on the surface properties, has been discussed in detail by Iler1 and Fubini.7 Biomolecule-silica interactions involve electrostatic, hydrogen bonding, and hydrophobic forces. The complexity and nature of these forces and their effect on equilibrium and nonequilibrium properties of biomolecules is poorly understood at the atomistic scale. A deeper understanding of the interactions between water, biomolecules, and amorphous silica using molecular simulations may provide further insight into physical, chemical, and toxic effects of amorphous silica. Theoretical modeling of the applications listed above requires tractable models for the hydrated, hydroxylated silica surface with a feasible pathway for inclusion of biomolecules. In this work, we have taken a step toward that goal by formulating a model for the undissociated silica surface using data from ab initio quantum chemical data. Since we anticipate large-scale equilibrium and nonequilibrium studies of biomolecules interacting with the silica surface, we have constructed the model as an extension of models that are well-calibrated for bulk silica and water. This is an advantage in terms of ensuring the desired structural and transport properties of the bulk region, and also facilitates interfacing our model with current force fields for biological molecules. However, the simplicity of our model and the fact that we retained the site charges of the BKS and SPC/E model constrain the flexibility of our model. This is evidenced in Figure 6 where a simultaneous description of the binding of water to a silanol group along various paths of approach was not achieved. We chose to match the most probable path of approach. If we would have weighted less favorable approaches more in our fit, it would have the effect of reducing our calculated heat of immersion. Because of the affinity of water molecules for silanol groups, we found that the local hydrophobic and hydrophilic character of the hydroxylated silica surface was connected to the local surface density of silanol groups. The simultaneous existence of hydrophobic and hydrophilic regions likely plays a role in determining the special adsorption characteristics of silica. For example, the adsorption of methanol is less affected by dehydroxylation than the adsoprtion of water, implying that hydrophobic interactions of siloxane bridges with methyl groups is an additional mechanism of adsorption.102 This effect could be analyzed using the interaction model developed in this work. There are several important directions needed for future refinement and extension of our model. It can be used to generate starting configurations for ab initio simulations of hydrated, hydroxylated silica. Studies of this type can reveal shortcomings in the local structure, e.g., bond lengths and angles, predicted by a model like ours. However, in our model, important structural features like the pattern of hydroxylation are imposed on the basis of what is known about the hydroxylated surface. Testing of our hydroxylation protocol by comparison with ab initio simulation or experiment would be very useful, but unfortunately, the time and system size required for a meaningful ab initio annealing and sampling of the surface may not be currently feasible. Experiments do furnish important guidance, but there is considerable variation in the surface structure implied by different experiments. Another important extension which is currently under development in our group is the inclusion of dissociated silanol groups. This is essential for the description of the surface charge density of the hydrated silica surface and the modeling of electrokinetic phenomena in this system. Acknowledgment. This material is based upon work supported by the National Science Foundation under Grant No.

11192 J. Phys. Chem. B, Vol. 111, No. 38, 2007 EEC-0425626. The calculations reported here were made possible by a grant of resources from the Ohio Supercomputer Center. Yun Kyung Shin is gratefully acknowledged for helpful discussions and assistance with hydroxylation and annealing protocol. References and Notes (1) Iler, R. K. The Chemistry of Silica; Wiley: New York, 1979. (2) Legrand, A. P., Ed. The Surface Properties of Silicas; Wiley: New York, 1998. (3) Helms, C. R., Deal, B. E., Eds. The Physics and Chemistry of SiO2 and the Si-SiO2 Interface; Plenum: New York, 1988. (4) Chabal, Y. J., Ed. Fundamental Aspects of Silicon Oxidation; Springer: New York, 2001. (5) Nawrocki, J. J. Chromatogr. 1997, A779 (1-2), 29-71. (6) Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; Mackerell, A. D., Jr. J. Phys. Chem. B 2006, 110 (6), 2782. (7) Fubini, B. In The Surface Properties of Silicas; Legrand, A. P., Ed.; Wiley: New York, 1998; p 415. (8) Woodcock, L. V.; Angell, C. A.; Cheeseman, P. J. Chem. Phys. 1976, 65 (4), 1565. (9) Tsuneyuki, S.; Tsukada, M.; Aoki, H.; Matsui, Y. Phys. ReV. Lett. 1988, 61 (7), 869. (10) Vashishta, P.; Kalia, R. K.; Rino, J. P.; Ebbsjo¨, I. Phys. ReV. 1990, B41 (17), 12197. (11) Ogawa, H.; Shiraishi, Y.; Kawamura, K.; Yokokawa, T. J. NonCryst. Solids 1990, 119 (2), 151. (12) Okuno, M.; Kawamura, K. J. Non-Cryst. Solids 1995, 191 (3), 249. (13) Feuston, B. P.; Garofalini, S. H. J. Chem. Phys. 1988, 89 (9), 5818. (14) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Phys. ReV. Lett. 1990, 64 (16), 1955-1958. (15) Watanabe, T.; Fujiwara, H.; Noguchi, H.; Hoshino, T.; Ohdomari, I. Jpn. J. Appl. Phys. 1999, Pt.2, 38 (4A), L366. (16) Cagin, T.; Demiralp, E.; Goddard, W. A., III. Phys. ReV. Lett. 1999, 82 (8), 1708. (17) Flikkema, E., S. B. Chem. Phys. Lett. 2003, 378 (5-6), 622. (18) Sarnthein, J.; Pasquarello, A.; Car, R. Phys. ReV. Lett. 1995, 74 (23), 4682. (19) Sarnthein, J.; Pasquarello, A.; Car, R. Phys. ReV. 1995, B52 (17), 12690. (20) Pasquarello, A.; Car, R. Phys. ReV. Lett. 1997, 79 (9), 1766. (21) Pasquarello, A.; Car, R. Phys. ReV. Lett. 1998, 80 (23), 5145. (22) Pasquarello, A.; Sarnthein, J.; Car, R. Phys. ReV. 1998, B57 (22), 14133. (23) Trave, A.; Tangney, P.; Scandolo, S.; Pasquarello, A.; Car, R. Phys. ReV. Lett. 2002, 89 (24), 245504. (24) Benoit, M.; Ispas, S.; Tuckerman, M. E. Phys. ReV. 2001, B64 (22), 224205. (25) Rahmani, A.; Benoit, M.; Benoit, C. Phys. ReV. 2003, B68 (18), 184202. (26) Van Ginhoven, R. M.; Jo´nsson, H.; Corrales, L. R. Phys. ReV. 2005, B71, 24208. (27) Huff, N. T.; Demiralp, E.; Cagin, T.; Goddard, W. A., III. J. NonCryst. Solids 1999, 253 (1-3), 133. (28) Ceresoli, D.; Bernasconi, M.; Iarlori, S.; Parrinello, M.; Tosatti, E. Phys. ReV. Lett. 2000, 84 (17), 3887-3890. (29) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91 (24), 6269-6271. (30) Feuston, B. P.; Garofalini, S. H. J. Chem. Phys. 1989, 91 (1), 564. (31) Feuston, B. P.; Garofalini, S. H. J. Appl. Phys. 1990, 68 (9), 4830. (32) Feuston, B. P.; Garofalini, S. H. J. Phys. Chem. 1990, 94 (13), 5351. (33) Litton, D. A.; Garofalini, S. H. J. Appl. Phys. 2001, 89 (11), 6013. (34) Rustad, J. R.; Hay, B. P. Geochim. Cosmochim. Acta 1995, 59 (7), 1251-1257. (35) Rustad, J. R.; Wasserman, E.; Felmy, A. R.; Wilke, C. J. Colloid. Interface Sci. 1998, 198 (1), 119-129. (36) Bakaev, V. A.; Steele, W. A. J. Chem. Phys. 1999, 111 (21), 9803. (37) Leed, E. A.; Pantano, C. G. J. Non-Cryst. Solids 2003, 325 (1-3), 48. (38) Cruz-Chu, E. R.; Aksimentiev, A.; Schulten, K. J. Phys. Chem. B 2006, 110 (43), 21497-21508. (39) Vessal, B. J. Non-Cryst. Solids 1994, 177 (1), 103-124. (40) Smith, W.; Greaves, G. N.; Gillan, M. J. J. Chem. Phys. 1995, 103 (8), 3091-3097. (41) Guillot, B.; Guissani, Y. Mol. Simul. 1997, 20 (1-2), 41. (42) Møller, C.; Plesset, M. Phys. ReV. 1934, 46 (7), 618-622. (43) Bartlett, R. J.; Silver, D. M. Int. J. Quantum Chem. Symp. 1974, 8, 271-276.

Hassanali and Singer (44) Binkley, J. S.; Pople, J. A. Int. J. Quantum Chem. 1975, 9 (2), 229-236. (45) Bartlett, R. J.; Silver, D. M. Int. J. Quantum Chem. Symp. 1975, 9, 183-198. (46) Pople, J. A.; Binkley, J. S.; Seeger, R. Int. J. Quantum Chem. Symp. 1976, 10, 1-19. (47) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (48) Ercolessi, F.; Adams, J. B. Europhys. Lett. 1994, 26 (8), 583. (49) Tabacchi, G.; Mundy, C. J.; Hutter, J.; Parrinello, M. J. Chem. Phys. 2002, 117 (4), 1416. (50) Izvekhov, S.; Parinello, M.; Burnham, J. C.; Voth, A. G. J. Chem. Phys. 2004, 120 (23), 10986. (51) Ma, Y.; Foster, A. S.; Nieminen, R. M. J. Chem. Phys. 2005, 122 (14), 144709. (52) Saengsawang, O.; Remsungnen, T.; Frizsche, S.; Haberlandt, R.; Hannongbua, S. J. Phys. Chem. B 2005, 109 (12), 5684-5690. (53) Ferrari, A. M.; Ugliengo, P.; Garrone, E. J. Phys. Chem 1993, 97 (11), 2671. (54) Sauer, J. J. Phys. Chem. 1987, 91 (9), 2315-2319. (55) Sauer, J.; Ugliengo, P.; Garrone, E.; Saunders, V. R. Chem. ReV. 1994, 94 (7), 2095-2160. (56) Smith, W.; Forester, T. R. J. Mol. Graphics 1996, 14 (3), 136141. (57) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98 (12), 10089-10092. (58) Perera, U. E. L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103 (19), 8577-8593. (59) Bakos, T.; Rashkeev, S. N.; Pantiledes, S. T. Phys. ReV. Lett. 2002, 88 (5), 55508. (60) Masini, P.; Bernasconi, M. J. Phys.: Condens. Matter 2002, 14 (16), 4133-4144. (61) Du, M.-H.; Kolchin, A.; Cheng, H.-P. J. Chem. Phys. 2003, 119 (13), 6418-6422. (62) Mischler, C.; Horbach, J.; Kob, W.; Binder, K. J. Phys.: Condens. Matter 2005, 17 (26), 4005-4013. (63) For example, see p 626 of ref 1. (64) Page 662 of ref 1, attributed to Kirichenko and Vysotskii.103 (65) Chuang, I.-S.; Kinney, D. R.; Maciel, G. E. J. Am. Chem. Soc. 1993, 115 (19), 8695-8705. (66) Chuang, I.-S.; Maciel, G. E. J. Am. Chem. Soc. 1996, 118 (2), 401406. (67) Chuang, I.-S.; Maciel, G. E. J. Phys. Chem. B 1997, 101 (16), 3052-3064. (68) Ong, S.; Zhao, X.; Eisenthal, K. B. Chem. Phys. Lett. 1992, 191 (3-4), 327-335. (69) Fisk, J. D.; Batten, R.; Jones, G.; O’Reilly, J. P.; M. Shaw, A. J. Phys. Chem. B 2005, 109 (30), 14475-14480. (70) Duval, Y.; Mielczarski, J. A.; Pokrovsky, O. S.; Mielczarski, E.; Ehrhardt, J. J. J. Phys. Chem. B 2002, 106 (11), 2937. (71) Rignanese, G.-M.; Charlier, J.-C.; Gonze, X. Phys. Chem. Chem. Phys. 2004, 6 (8), 1920-1925. (72) Walsh, T. R.; Wilson, M.; Sutton, A. P. J. Chem. Phys. 2000, 113 (20), 9191-9201. (73) Murashov, V. V. J. Phys. Chem. B 2005, 109 (9), 4144. (74) Zhuravlev, L. T. Pure Appl. Chem. 1989, 61 (11), 1969-1976. (75) Rarivomanantsoa, M.; Jund, P.; Jullien, R. J. Phys.: Condens. Matter 2001, 13 (31), 6707-6718. (76) Roder, A.; Kob, W.; Binder, K. J. Chem. Phys. 2001, 114 (17), 7602-7614. (77) Taylor, J. A. G.; Hockey, J. A. J. Phys. Chem. 1966, 70 (7), 21692172. (78) Taylor, J. A. G.; Hockey, J. A.; Pethica, B. A. Proc. Brit. Ceram. Soc. 1965, 5, 133-141. (79) Makrides, A. C.; Hackerman, N. J. Phys. Chem. 1959, 63 (4), 594598. (80) Takei, T.; M. Chikazawa. J. Colloid Interface Sci. 1998, 208 (2), 570-574.

Model for the Water-Amorphous Silica Interface (81) Takei, T.; Eriguchi, E.; Fuji, M.; Watanabe, T.; Chikazawa, M. Themochim. Acta 1998, 308 (1-2), 139-145. (82) Flyvbjerg, H.; Petersen, H. G. J. Chem. Phys. 1989, 91 (1), 461. (83) Warne, M. R.; Allan, N. L.; Cosgrove, T. Phys. Chem. Chem. Phys. 2000, 2 (10), 3663. (84) Fubini, B.; Bolis, V.; Cavengo, A.; Garrone, E.; Ugliengo, P. Langmuir 1993, 9, 2712. (85) Bolis, V.; Cavenago, A.; Fubini, B. Langmuir 1997, 13 (5), 895. (86) Israelachvili, J. Intermolecular and surface forces, 2nd ed.; Academic: San Diego, CA, 1991. (87) Lyklema, J. Fundamentals of Interface and Colloid Science; Academic: San Diego, CA, 1991. (88) Jal, P.; Patel, S.; Mishra, B. Talanta 2004, 62 (5), 1005. (89) Conlisk, A. T.; Singer, S. J. In Biomolecular Sensing, Processing and Analysis; Bashir, R., Wereley, S., Ferrari, M., Eds.; Handbook of BIOMEMS and Biomedical Nanotechnology, Vol. 4; Springer: New York, 2006; p 301. (90) Cordek, J.; Wang, X.; Tan, W. Anal. Chem. 1999, 71 (8), 1529. (91) Liu, X.; Tan, W. Anal. Chem. 1999, 71 (22), 5054. (92) Fang, X.; Liu, X.; Schuster, S.; Tan, W. J. Am. Chem. Soc. 1999, 121 (12), 2921.

J. Phys. Chem. B, Vol. 111, No. 38, 2007 11193 (93) Ingersoll, C. M.; Jordan, J. D.; Bright, F. V. Anal. Chem. 1996, 68 (18), 3194. (94) Hartnett, A. M.; Ingersoll, C. M.; Baker, G. A.; Bright, F. V. Anal. Chem. 1999, 71 (6), 1215. (95) Wang, A. J.; Arnold, M. A. Anal. Chem. 1992, 64 (9), 1051. (96) Kar, S.; Arnold, M. A. Anal. Chem. 1992, 64 (20), 2438. (97) Valle-Delgado, J. J.; Molina-Bolı´var, J. A.; Galisteo-Gonza´lez, F.; Ga´lvez-Ruiz, M. J.; Feiler, A.; Rutland, M. W. Langmuir 2005, 21 (21), 9544. (98) Lundqvist, M.; Andresen, C.; Christensson, S.; Johansson, S.; Karlsson, M.; Broo, K.; Jonsson, B.-H. Langmuir 2005, 21 (25), 11903. (99) van der Veen, M.; Norde, W.; Stuart, M. C. Colloids Surf. 2004, B35 (1), 33. (100) Sui, J.; Tleugabulova, D.; Brennan, J. D. Langmuir 2005, 21 (11), 4996. (101) Hassanali, A. A.; Li, T.; Zhong, D.; Singer, S. J. J. Phys. Chem. B 2006, 110 (21), 10497. (102) Page 653 of ref 1. (103) Kirichenko, L.; Vysotskii, Z. Z. Dokl. Akad. Nauk SSSR 1967, 175 (3), 635.

Model for the WaterâAmorphous Silica Interface: The Undissociated

Recommend Documents

Model for the WaterâAmorphous Silica Interface: The Undissociated