Choice of Computational Techniques and Molecular Models for ab

Daniel Bougeard, Konstantin S. Smirnov, and Ekkehard Geidel ... Brian J. Teppen, Kjeld Rasmussen, Paul M. Bertsch, David M. Miller, and Lothar Schäfe...
0 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. 1994, 98, 12545-12557

12545

Choice of Computational Techniques and Molecular Models for ab Znitio Calculations Pertaining to Solid Silicates Brian J. Teppen,? David M. Miller,? Susan Q. Newton: and Lothar Schdifer*,* Department of Agronomy and Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701 Received: March 14, 1994; In Final Form: July 27, 1994@

Development of potential energy parameters for the modeling of condensed-phase silicates is relying increasingly on the results of a b initio calculations. We examined the effects of basis set size and electron correlation corrections on the properties of three increasingly complex molecules often used as models for and disilicic acid (H6Si207). When disiloxane silicates, disiloxane (H3SiOSiH3), orthosilicic acid was optimized with the MC6-311G(2d,2p) basis, MP2 corrections changed the HF results through narrowing the Si-0-Si angle by 8" and lengthening each Si-0 bond by 0.022 A. When the same basis was applied to Si(OH)4, the MP2-optimized bond lengths for Si-0 and 0 - H were 0.020 and 0.017 A larger, respectively, than their HF-optimized values. Meanwhile, the Si-0-H angle decreased 3" from its HF value, to 115.4' at the MP2 level. In comparing these results to the HF/and MP2/MC6-311G** optima, it was evident that the basis set limit had not been reached. Due to the hydrogen-bonding interactions in H6Si207, the MP2/ MC6-3 l l G * * correlation corrections to geometries of the three conformers considered here were more dramatic. The Si-0 bond lengths increased by as much as 0.029 A over their HF/MC6-311G**-optimized values, 0 - H bonds lengthened by up to 0.021 A, and Si-0-Si angles decreased by up to 13.7'. Some framework Si-0-Si-0 torsions changed by 10" when electron correlation effects were included in the optimization, and nonbonded OH distances decreased by as much as 0.62 A. Further optimizations showed that only two of the conformers were stable at the MP2 level, emphasizing the importance of electron correlation in the study of silicate molecular models. The two MPZoptimized structures of H6Si207 were separated by 2.4 kcal mol-' , but when zero-point energies and thermal corrections were included, the conformers were within 0.4 kcal mol-'. It is argued on the basis of electrostatic potential maps and a b initio partial charges that H3SiOSiH3 and its analogues provide poor quantitative models for silicate functional groups. However, structural features and partial charges on Si(OH)4 and H6Si207 proved to be very similar at a given level of theory. Thus, effective a b initio parameter development for silicates should focus on molecules comprised of silicon that is fully coordinated by oxygen.

Introduction

with the crystal lattice. However, the successes of the molecular approach have been attributed to the overriding importance of There is tremendous interest in modeling sorption and local geometry trends in establishing crystal structure^.^^^^^-^^ catalysis processes in zeolites and clays, in order to improve Indeed, partial charges and geometries from ab initio calculaour understanding of natural processes and current technologies tions on the H6Si207 molecular fragmentz1agreed rather well and also to design better sorbents and catalysts. The most with those calculated using the same minimal basis set and the important prerequisite to simulating such materials by molecular full translational symmetry of either the quartz crysta121 or dynamics or Monte Carlo methods is the development of siliceous zeolite^.^^^^^ The optimized geometries of molecular accurate interatomic potential energy functions. While many clusters seem to be retained in solid silicates because the highly parameters are needed to objectify these functions, relatively flexible Si-0-Si angle bends to accommodate strain with very few experimental data are applicable, so silicate researchers are little energetic penalty.z8~38~39 forced to follow Clementi's example' and gather conformational Since local geometry trends have been shown to be so crucial potential energy data through appropriate ab initio cal~dations.~-'~ Solid silicates are essentially infinite systems when viewed for predicting silicate crystal structures, one would expect that from the atomic perspective. Truly ab initio methods exist18-z0 ab initio geometry optimizations including the effects of electron to solve the full crystal wave function using the translational correlation would be necessary in order to develop accurate symmetry of the solid, but for complex silicates, only minimal potential energy parameter~.~O-~* Ab initio based silicate basis sets may be used and even then only a few degrees of parameter development that is not grounded in calculations freedom optimized even on supercomputers.z1-26 Therefore, which include electron correlation cannot arrive at self-consistent ab initio calculations pertaining to solid silicates have focused parameters for electrostatic or dispersion interactions." The on small silicate molecules representative of crystal substrucobjective of this research, therefore, was to systematically tures; such calculations have been reviewed several time^.^'-^^ examine the effects of electron correlation on the predicted Caution is warranted in using this approach,34 because geometries of three small silicate molecules using two robust, calculations on isolated molecules ignore cumulative interactions polarized basis sets. The molecules we chose to study, H3SiOSiH3, Si(OH)4, and (H0)3SiOSi(OH)3, have often been Department of Agronomy. used as models for silicate functional groups and exhibit an * Department of Chemistry and Biochemistry. orderly increase in size and complexity. This hierarchic * Abstract published in Advance ACS Abstracts, November 1, 1994. 0022-365419412098-12545$04.50/0

0 1994 American Chemical Society

12546 J. Phys. Chem., Vol. 98, No. 48, 1994 a p p r ~ a c hallows ~ ~ . ~comparison ~ of the relative suitabilities of these molecules as models for solid silicates.

Teppen et al.

tion functions on all atoms. Very similar results, including a linear structure for disiloxane, were found by Nicholas and coworkers,12 who used a 'TZ+d" basis set that differed from Methods MC6-311G** only by using a Dunning sp basis7I on 0 (1%3pl) and H ([3s]) rather than 6-311G. All ab initio geometry optimizations were performed using When the HF/MC6-311G** results were improved through the Gaussian92 system of programs4 to calculate single-point the addition of MP2 perturbation corrections for electron wave functions, energies, population analyses, and forces, while correlation, the Si-0 bonds in disiloxane lengthened by 2.0 the forces were relaxed in normal coordinate space45according pm as one might expect. The greatest effect was seen in the to the B-matrix gradient method of Pulay and c o - w ~ r k e r s . ~ , ~ ~ Si-0-Si angle, which closed from 180" to 157.6'. The 22' This procedure is similar to that described by C6rsky et a1*! change is huge, but it is well-known that the potential energy and implemented by Sauefl3 for silicate calculations. All surface for Si-0-Si bending in disiloxane is very flat, with geometry optimizations were completely unconstrained (that is, an angle inversion barrier estimated to be only about 112 cm-1?7 only C1 symmetry was assumed). An optimized geometry was just 0.32 kcal mol-'. Similarly, large changes in the Si-0-Si defined as one for which all internal coordinate forces were angle when electron correlation is included have been seen with less than 0.001 mdyn. The calculations were done by an IBM other basis sets (see Table 1). RSf6000 Model 350 workstation. Basis sets used on 0 and H atoms were 6-311G**,49 but Based on this evidence alone, the inclusion of electron the basis for Si was derived from the research of McLean and correlation would seem mandatory to adequately describe Si-0 Chandler.5o The latter basis is contracted as [6s,5p]: The (12s,bonding. However, an unsettling fact about this situation is 9p) primitives are distributed as (631 111,4211l), giving a basis that the 6-31G basis,56 although smaller than MC6-311G, set for Si that is fully double-t and is triple-f in the 3p orbitals. predicts a bent disiloxane geometry at the HF level when This hybrid basis set has been recommended5' as a starting point augmented with polarization function^^^,^^ on just heavy atfor large-basis-set correlated calculations on molecules containo m ~or ~on all ~ ' ~ ~ It is perhaps noteworthy that the ing second-row elements. Polarization functions were uncon6-31G* basis uses a set of six d functions to allow for tracted and were deployed in two different manners. In the fust polarization, rather than the five "pure" d functions included in basis set, hereafter denoted MC6-311G**, one set of five d the larger basis sets of K0put,6~Nicholas et al.,12 and the present functions was used on each Si or 0 atom, with d exponents of study. The extra d function may allow the wave function more 0.450 on Si52and 1.292 on 0,49 while the p exponent on H was flexibility about the Si and 0 nuclei. In a cryptic note in one 0.750!9 In the second basis set, which we label MC6-311Gof their papers, Gibbs and c o - ~ o r k e r sindicate ~~ that if three (2d,2p), double-f polarization functions were found to be useful; sets of d functions are used in conjunction with the 6-31G the two sets of five d functions for Si had e ~ e n - t e m p e r e d ~ ~ ~basis, ~ ~ then the Si-0-Si angle in disiloxane optimizes to 143.7", exponents of 0.900 and 0.225, the two d exponents for 0 were very close to the gas electron diffraction value@ of 144.1'. 2.584 and 0.646, and the two p exponents for H were 1.500 Ugliengo and c o - ~ o r k e r sdetermined ~ ~ ~ ~ ~ that the greatest and 0.375. For one optimization, we used the standard 6-31G* improvement in basis sets for silanol (SiH30H) came from basis set52156-57 in order to supplement literature calculations splitting the polarization functions into two independent sets. using 6-31G* that did not include electron correlation. Finally, the recent systematic investigation of Nicholas et al.l2 In all calculations that included correlation effects, we used showed that adding a second set of d functions on each heavy the second-order perturbation (MP2) method of Moller and atom allows a triple-f basis set to make adequate structural Plesset,s8 which has been shown to typically account for 85predictions about disiloxane at the HF level of quantum theory. 95% of the correlation energy in reasonable computational They also found that adding a third set of polarization functions time^.^^,^ Of course, variational methods such as full configugave little added benefit. ration interaction would be preferable to perturbative methods Another disturbing observation is that altering the value of but would be impossible for the systems considered here. In the d exponent can have large effects on the ab initio geometry addition, the size consistency60 of perturbative methods, as predictions for H3SiOSiH3. Indeed, even the small 3-21G* opposed to incomplete CI methods, is a desirable feature for basis set (which again uses six d functions) can predict a bent modeling purposes. Finally, MP perturbation theory is the disiloxane if the silicon d exponent is chosen to do In simplest treatment of dynamical electron correlation61 and, their study of Si and 0 polarization functions, Grigoras and therefore, the method most likely to be used by nonspecialists, Lane67 showed that it is possible to optimize the d exponent so the use of MP2 calculations will best allow us to compare for Si so that the molecular energy is minimized or so that some our results with those of other geochemists. All electrons were other properties of the molecule (Le., the Si-0-Si angle) are involved in the perturbation corrections. The MP2 population best reproduced. Earley70 showed that the Si-0-Si angle was analyses all included the true MP2 density62 based on the similarly sensitive to the value of the d exponent when the Z-~ector.~~ 6-31G* basis set is used. The Si d-orbital exponent (0.45) that we used in MC6-311G** is the one that minimizes the Results and Discussion energy of di~iloxane,~~ Si02,8a and an average of several Sicontaining molecules.52 However, Grigoras and Lane67found Disiloxane. The simplest neutral molecule containing an 0 that d exponents of 0.30 or 0.90 led to the most faithful bridging between two Si atoms is disiloxane, H3SiOSiH3. reproduction of the disiloxane geometry with the 3-21G* basis Experimental information indicates that disiloxane has a bent set. Their goal was a small basis set, so they used 0.3 as their framework geometry with an Si-0-Si angle of 140 to 150".64*65 best value. We had the luxury of adding a second set of Unfortunately, researcher^^^@'-^^ have found that the Si-0polarization functions on each atom, as recommended by Si angle in disiloxane optimizes to 180" With many different Nicholas and associates.l 2 Common practices5 when adding basis sets at the Hartree-Fock (HF) level of q antum theory. another set of polarization functions is to make the exponents As outlined in Table 1, this was also the c se for our HF follow a geometric progression separated by a factor of 4, which calculations using the MC6-311G** basis set, even though the means that 0.45 gives rise to 0.225 and 0.90. Thus, the lowest basis is of nearly valence-triple-f quality and includes polariza-

gd

J. Phys. Chem., Vol. 98, No. 48, 1994 12547

Silicate Molecular Models

TABLE 1: Results from the Present and Previous ub Initio Studies of Selected Properties for Disiloxane, as Well as Some Experimental Measurements of Those Properties methodbasis setlstudy MP2/MC6-311G(2d,2p) (present paper) MP2(frozen core)/6-31l G ( 2 d , ~ ) ~ ~ MP2/MC6-311G** (present paper) MP2(frozen core)/6-31lG**72 ClSDISi:(12s9p2d)/[6s5p2d] l2 0:(10s6p2d)/[5s3p2d] H:(%2p)/[3s2p] CPP3/Si:(11s7p 1d)/[6s4p ld]@ 0:(9~5pld)/[5~3pld] H:(5s lp)/[3s 1p] Mp2/6-3 1G*70.72,74 HF/Si:( 12~9p2d)/[6~5p2d]'~ 0:(10s6p2d)/[5s3p2d] HF/MC6-311G(2d,2p) (present paper) HF/MC6-311G** (present paper) HF/6-31G(3d)37 HF/Si:(1 1s7p 1d)/[6s4p ld]69 0:(9s5pld)/[5s3pld] H:(Sslp)/[3slp] HF/6-31G**75 ~ ~ / 61~*72,76 - 3 m ~ c p - 13~ * 7 0 HF/3-21G*68 HF/3-2 1G(*)o,67

Si-0, A

energy, hartrees Methods Including Electron Correlation -657.0772 -656.7807 -657.0052 -656.7378 -656.77 19

Si-0-Si, deg

1.645 1.647 1.640 1.641 1.630

139.3 139.8 157.6 156.5 143.8

-656.7330

1.649

150.4

-656.5916

1.656-1.660

143.6- 145.2

Methods Neglecting Electron Correlation -656.3558 1.620 -656.3467 -656.3344 -656.2854 -656.2819

-656.2590

147.3

1.623 1.621 1.632 1.627

147.1 180.0 143.7 180.0

1.625 1.626 1.646 1.630 1.645

160.0 162-174 143.4 180.0 149.5

1.634 f 0.002 1.631 f 0.006

144.1 f 0.9 142.2 f 0.3

Experimental Methods gas electron diffractionM low-temp X-ray d i f f r a ~ t i o n ~ ~ (I

The d exponent on Si (0.30) was optimized to reproduce the experimental Si-0-Si

TABLE 2: Optimized Geometriep Derived from the ab Initio Study of Disiloxane propertyb

HF/MC6-311G(2d,2p) MP2/MC6-31lG(2d,2p) 1.623 1.623 1.469 1.473 1.473 1.473 1.473 1.469 147.1 108.7 110.1 110.1 110.1 110.1 108.7 -60.3 59.7 179.7 -179.7 60.3 -59.7

1.645 1.645 1.468 1.473 1.473 1.473 1.473 1.468 139.3 108.1 110.3 110.3 110.3 110.3 108.1 -60.4 59.8 179.7 -179.7 60.4 -59.8

Bond lengths in angstrams; bond angles and torsions in degrees.

* Atom numberings taken from Figure 1.

energy d exponent52splits to produce approximately the two d exponents that produce the best Si-0-Si g e ~ m e t r i e s .When ~~ this new basis set (MC6-311G(2d,2p)) was applied to disiloxane, the results were striking, as should have been expected after the work of Nicholas et al.12 The HF/MC6-311G(2d,2p) energy minimized at an Si-0-Si angle of 147.1°, rather than the linear structure found with MC6-311G**. This result, along with other geometric details presented in Tables 1 and 2, is in good agreement with the calculations of Nicholas et al.12 using their "TZf2d" basis and is in reasonable agreement with the gas electron diffraction data.64

angle.

Since increasing the number of polarization functions in the basis set allows the Hartree-Fock results to improve so much, it is of interest to investigate whether electron correlation corrections remain significant for disiloxane. A full MP2/ MC6-3 11G(2d,2p) optimization was performed; the resultant MP2 electron density and its difference from the HF density (at the MP2 geometry) are plotted in Figure 2. The inclusion of electron correlation allowed a significant increase in the electron density immediately adjacent to each nucleus, with larger increases near Si than 0. Also, electron, density was removed from the higher-density bonding regions but added to the lower-density bonding regions and to the very diffuse regions of the molecule. In this sense, electron correlation decreased the ionic character of the Si-0 bond by reinforcing the weakest part of that bond. Geometric effects of correlation (Table 2) included the lengthening of Si-0 bonds by a rather substantial 2.2 pm, while the Si-H bonds contracted slightly. These changes are similar to those observed in correcting MC6-3 llG** by MP2, as well as aose found by K ~ p uwhen t ~ ~CPP3 corrections were applied, except that the Si-H bond lengthened slightly in his study. Curtiss et al.,74Earley,'O and Luke72 found an even larger increase of 3.0-3.2 pm in the Si-0 length when MP2 optimizations were done using the 6-31G* basis. This agrees with observationsmsg1 that the less flexible the basis set, the more overexaggerated will be the MP2 corrections to the geometry. Nicholas et al.'* found a much smaller increase, 1.0 pm, in the Si-0 bond length when correlation effects were included at the CISD level of theory. Their result for the correlated bond length agrees more closely with the experimental determinat i o n ~ than ~ , ~do~our calculations. Those same researchers12 found that CISD corrections decreased the Si-0-Si angle by 3.5", while our MP2 results show a decrease of 7.6" with a similar basis set. Again, this substantial change brings our MP2

Teppen et al.

12548 J. Phys. Chem., Vol. 98, No. 48, 1994

disiloxane disilicic acid, conformer A

silicic acid

disilicic acid, conformer B

W

disilicic acid, conformer C

Figure 1. Atom numbering for the energy minima of disiloxane, orthosilicic acid, and three conformers of disilicic acid.

calculations into poorer agreement with experiment64.65than the CISD calculation. Poorer agreement with some experimental data does not, of course, invalidate our MP2 results. Nicholas and associates12 rightly point out that ab initio results for the very flexible Si0-Si angle should not be compared directly with experiment because thermal excitation of higher vibrational states contributes to the measured angle. They estimate that Boltzmann weighting would increase the room-temperature Si-0-Si by about 5.5' above its ground-state value, thereby indicating that our MP2/MC6-31 1G(2d,2p) result of 139.3' would predict an experimental angle of 144.8', in very good agreement with the gas electron diffraction data.a This is also consistent with the finding65thaf the Si-0-Si angle in crystalline disiloxane has decreased to 142.2" at 108 K, although the crystalline environment may distort the angle somewhat. The well-known work of Gibbs27*28*37.82-84 has shown that there is a nonlinear inverse relationship between the Si-0-Si angle and the Si-0 bond length. Thus, the increase in Si-0-Si due to thermal averaging would cause the experimental Si-0 bond length to decrease compared to the ground state, possibly explaining why our M p u MC6-311G(2d92p) bond length is 0.011 A longer than the experimentala value. By the same token, the variational theorem does not apply to MP2 calculations, so the substantially lower energy of our structure does not necessarily imply a better geometry. However, the trend in the systematic study of Nicholas et al.12 was for each improvement in the quality of the calculation to decrease Si-0-Si: Adding a second set of polarization functions and increasing the number of s and p functions and CISD correlation corrections all moved the equilibrium angle to smaller values. Our MP2 calculations decrease that value even farther. Furthermore, the CISD-optimized structure12was an extrapolation along the potential energy surface probed by a series of single-point calculations,rather than a full optimization, and core electrons were not allowed to contribute to the excited states. Our MP2 calculations, on the other hand, included all electrons and were unconstrained geometry optimizations. It

Figure 2. Total MP2/MC6-311G(2d92p) electron density for H3SiOSiH3 (top) and density difference map in the Si-0-Si plane (bottom) at the MF'2-optimizedgeometry for I-hSiOSiHJ. The difference map illustrates the effects of electron correlation on the electron density by subtracting the total HF density from the total MP2 density. The numerical density values are given in atomic units.

is typical62of correlation correction studies that MP2 optimizations somewhat overexaggerate correlation effects while CISD optimizations underestimate them. It would be of interest to extend such an optimization to the MP4/MC6-31 IG(2d-2~) level of perturbation theory to assess whether the MP2 corrections in this case have gone too far or not far enough. Orthosilicic acid. The basic structural unit of most silicate minerals is the Si tetrahedron, and the simplest molecule containing this tetrahedron is silicic acid, Si(OH)4. The point group symmetry for the global energy minimum of silicic acid is The initial geometries used in our optimizations were nearly S4, but no symmetry constraints were enforced. During all optimizations, approximate S4 geometries were retained. Figure 1 shows the molecular geometry and Table 3 presents comparisons of energies and geometries for Si(OH)4 when calculated with and without the inclusion of electron correlation effects. Of the four basis sets that permit direct comparison of HF and MP2 results, namely, MC6-311G(2d92p) and MC6-311G** (present work), the Dunning (9sSp) basis,87 and 6-31G*43,86*88(present work for M E ) , the first two are the largest and consequently give the lowest HF energies for Si(OH)4. Note that, even with rather robust basis sets, both S4?96737.439&1-87

J. Phys. Chem., Vol. 98, No. 48, 1994 12549

Silicate Molecular Models TABLE 3: Selected Properties Derived from the ab Initio Study of the with Some Previous Results

methodhasis set/study MP2/MC6-311G(2d,2p)(present paper) MP2/MC6-311G** (present paper) MP2lDunning (9~,5p)**~’ MP2/6-31G* (present paper) HF/MC6-311G(2d,2p)(present paper) HF/MC6-311G** (presentpaper) HFlDunning (9~,5p)**~~ HF/6-31G**6 HF/6-31G*43.86.88

Sq

Conformation of Orthosilicic Acid and Comparison

Si-0, A 1.634 1.641 1.658 1.651 1.614 1.622 1.632 1.621 1.629

energy, hartrees -592.2471 -592.1 172 -591.9282 -591.7153 -591.0557 -591.0243 -590.9973 -590.89

0-H, 8, 0.955 0.957 0.964 0.970 0.938 0.939 0.943 0.944 0.947

Si-0-H, deg 115.4 114.7 114.7 114.2 118.4 119.5 118.5 130.1 117.1

TABLE 4: Optimized Geometriefl Derived from the ab Initio Study of the S4 Conformation of Si(OH)4

propertyb Si-0 (all) 0-H (all)

02-Si[-03 02-Sil-04 02-Sil-05 03-Sil-04 03-Si[-05 04-Sil-05 Sil-02-& Sil-03-H7 Sil-04-H8 Si1-05-H9 OZ-Sil-03-H7 03-Sil-05-Hg 04-Sii Os-Sil-04-H8 a

HF/MC6-311G** 1.622 0.939 106.8 106.7 115.1 115.1 106.7 106.7 119.4 119.5 119.5 119.5 36.2 -36.2 -35.9 36.2

MPuMC6-311G** 1.641 0.957 106.2 106.2 116.3 116.3 106.2 106.2 114.7 114.7 114.7 114.7 37.0 -36.9 -37.2 36.8

HF/MC6-311G(2d,2p) 1.614 0.938 106.8 106.7 115.1 115.1 106.7 106.7 118.3 118.4 118.4 118.4 36.8 -36.8 -36.4 36.6

MP2/MC6-31lG(2d,2p) 1.634 0.955 106.4 106.3 115.9 115.8 106.4 106.4 115.4 115.4 115.4 115.4 37.2 -37.2 -37.2 37.0

Bond lengths in angstroms; bond angles and torsions in degrees. Atom numberings taken from Figure 1.

the HF and the MP2 results remain basis set dependent, with Si-0 bond lengths, for example, ranging from 1.614 to 1.632 A at the HF level and 1.634 to 1.658 8, at the MP2 level. The shortest bond lengths in each set arise from the largest calculations. Since Nicholas et al.12 have already shown the importance of a second set of polarization functions for H3SiOSiH3, Table 4 displays the MC6-311G** and MC6-311G(2d,2p) geometries of Si(OH)4 for detailed comparison. At the HartreeFock level, extra polarization functions normally decrease bond lengths by decreasing two-electron repulsions and thus allowing an increase in the wave function amplitude between the n u ~ l e i This . ~ ~behavior ~ ~ ~ is borne out for Si(OH)4, as the second set of polarization functions shrinks the Si-0 bonds by 0.8 pm and the 0-H bonds shrink slightly. Another effect of increased polarization on Hartree-Fock wave funtions, especially at electron-rich nuclei such as 0,81,89 is to allow better angular flexibility. Thus, the HF/MC6-311G(2d,2p) values for the Si-0-H angles move toward their correlated values as compared to the HF/MC6-311G** angles. De Almeida and O’Malley86~90-92 have repeatedly asserted that “the inclusion of electron correlation at the MP2 ...level of theory have (sic) little effect on the values of geometric parameter^"^^^^^ for silicic acid,s6 silicate ion,91and silan01.~~ The evidence they cite is that, for the silicate ion [H3Si04]- at the MP2/6-31G* level of theory, Si-0 bond lengths increased “only” 3.0 pm, 0-H bond lengths increased “only” 2.3 pm, and bond angles changed less than 4” as compared to the HF/ 6-3 1G*-minimized structure.91 They recommend MP2 singlepoint calculations at HF/6-3 1G*-minimized geometries as an adequate theoretical approachg1and comment that “the HF/631G* treatment is shown to be adequate for predicting the geometry of the [H3Si04]- silicate species.”91 Nevertheless, the greatest effect of electron correlation would be expected for nonbonded interactions, and torsional angles controlling intramolecular hydrogen bonding in [H3Si04]-

changed by as much as 11” when MP2 corrections were applied?l In addition, those same authors observed that electron correlation using MP2 techniques decreased the vibrational frequencies by 20-25%86 and increased internuclear distances of transition-state structures by When MP2 theory is applied to silicic acid using the four basis sets outlined in Table 3, Si-0 and 0-H bonds do increase by “only” 2-3 pm, which amounts to more than 2% of the 0-H bond length. The Si-0-H angles change by 3-5” and the 0-Si-0 angles by about 1”. Ferriss7has shown that, taken together, these effects have unexpected and nonsystematic consequences for the vibrational spectra of silicic acid. Both the order and relative intensities of the normal modes were different at the SCF and correlated levels of treatment, especially for the low frequenciess7 that have the greatest impact on molecular dynamics and thermodynamic predictions. Generally speaking, the sizes of the M p 2 geometry corrections (Table 3) are again smaller as the quality of the basis ‘set increases from 6-31G* to MC6-311G(2d,2p). Table 4 illustrates the individual and combined effects of adding more polarization functions and electron correlation corrections to Si(OH)4, the prototypical silicate tetrahedron. Since the MP2 wave function is dominated by the HF contribution, an extra set of polarization functions allows MP2-optimized bond lengths to decrease and angles to be more flexible. The added flexibility of the largest basis set means that the correlated Si-0-H angle of 115.4” should be more accurate than the 114.7” found by Ferriss7 and by the present study with less polarized basis sets. It is disturbing that the Si-0 bond length determined by the HF/MC6-3 11G(2d,2p) calculations agrees with the experimental bond length in a-quartz (1.60993to 1.610g4A) better than the results of any of the MP2 calculations. However, the MP2/ MC6-31 lG(2d,2p)-determined Si-0 distance of 1.634 8, is near the center of the range of Si-0 distances (1.55- 1.70 8,) found in tetrahedral silicate minerals as a who1e.s3~s4,95-97 Note that our highest-level MP2 optimizations of Si-0, 0-H, and

Teppen et al.

12550 J. Phys. Chem., Vol. 98, No. 48, 1994

Si-0-H agree very well with S a ~ e r ' sextrapolations ~~ from an HF/6-3 lG* calculation. Disilicic Acid. The simplest neutral molecule that allows us to study interactions between two Si tetrahedra is H6Si207, which contains the Si-0-Si linkage in the context of Si atoms that are fully coordinated by 0 atoms. The bridging 0 is now surrounded by two whole coordination shells as found in silicate minerals, so it should have a reasonably complete local bonding environment. Hydrogens are added to the terminal oxygens to better simulate the bonding environments of the other atoms.30 The addition of these H atoms, however, complicates the molecular potential energy surface by allowing for several pattems of hydrogen bonding to occur. After a number of semiempirical and STO-3G calculations on this hypothetical molecule, O'Keefe et al.36 seem to have been the first to use an extended basis set (6-31G) to examine H6Si207. However, they constrained their geometries in several ways, as did Lasaga and G i b b ~ The . ~ ~latter authors late+ lifted all constraints, using the 6-31G basis with a set of d functions added to 0 atoms only. The optimized Si-0-Si angle for the structure they reported6 was 137.6". This is already in reasonable agreement with the angles observed in d i ~ i l o x a n and e~~~~ in quart^,^^,^^ even though no polarization functions were added on Si. The reasonable angles in disilicic acid are probably due to intramolecular hydrogen bonding, since Hill and Sauer17 recently showed that, if such H bonding is prevented, then HF optimization with a fairly large basis set results in a linear Si0-Si structure. The optimized H6Si207 structure found by Lasaga and Gibbs6 resembled conformer A of Figure 1, featuring a single hydrogen bridge between hydroxyls in the two Si tetrahedra. Burkhard et al.99used the 3-21G* basis set to examine this conformer and also a second conformer (B of Figure 1) with two intertetrahedral H bridges. They found that the second H bond of conformer B was able to compress the Si-0-Si angle to 122.7', compared to 135.4' in conformer A, further illustrating the shallowness of that angle's potential energy surface. The authors99noted that an H bond of 1.94 8, formed in the sing1 bridged conformer (conformer A) and that H bonds of 2.12 formed in the doubly bridged conformer (conformer B). Burkhard et al.99found conformer B to have lower energy, by 1.93 kcal mol-'. Kubicki and SykeslOOstudied conformers B and C (Figure 1), with the same basis set. Conformer C features one H bond that is stronger than that of conformer A and also a second, very weak, H bridge. Their research" showed both structures to be true minima, with conformer C (Si-0-Si = 141.5') just 0.09 kcal mol-' lower in energy than conformer B (Si-0-Si = 128.3'). We examined these three structures with much larger basis sets than any previous study and, for the first time, investigated the magnitude of electron correlation corrections. The equilibrium energies are presented in Table 5 , and the structural data are detailed in Table 6. The ab initio energy relationships found for the HF/MC631 1G** optimizations support the previous 3-21G* findi n g ~ The . ~HF/MC6-311G** ~ ~ ~ ~ vibrational analysis indicates that conformers A and C are minima, while the Hessian of conformer B has one negative eigenvalue even after very precise convergence of the optimization. The three lowest-frequency normal modes occur at 22, 36, and 48 cm-' for conformer A, -14, 44,and 55 cm-' for B, and 25, 40, and 79 cm-' for C. When conformer B is optimized with the 3-21G** basis set, no eigenvalues of the Hessian are negative," and the lowest non-zero frequencies are 56, 81, and 108 cm-'. The 81-cm-'

1

TABLE 5: (a) Ab Initio Energy Relationships among Three Conformers (See Figure 1) of Disilicic Acid, €I&i2O7, and (b) Scaleda HF/MC6-311G** Zero-Point Energies (ZPE), Entropies (ST), and Relative Thermal Corrections (H- Ho), Free Energies (AG) for Three Conformers of Disilicic Acid at 300 K HF/MC6-311G** energy

MPUMC6-311G** energy

abs, hartrees

rel, kcal mol-'

abs, hartrees

rel, kcal mol-'

1.37 0.08

-1107.958 97 -1 107.962 85

2.44

Cb

-1106.008 44 -1 106.010 49 -1106.01062

conformer

ZPE, kcal mol-'

kcal mol-I

62.354 62.822 62.816

70.409 70.609 70.645

conformer A

B

A Bd

Cb

0

0 H

- Ha,

ST,

AG:

kcal molw1 kcal mol-' 33.925 32.552 33.119

0.40

0

Thermodynamic functions were calculated using a 0.89 scaling factor for all ab initio vibrational frequencies. Conformer C was a true minimum at the HF level but collapsed to conformer B at the MP2 level. Thus, no AG was calculated for conformer C. The MP2NC6311G** energies were used to calculate these relative free energies. dConformer B possessed one imaginary frequency at the HFMC6311G** level. This normal mode is a wag of the bridging 0 and corresponds to a vibration at 81 cm-l when conformer B is optimized with the 3-21G** basis set. The ZPE, H - Ho, and ST for conformer B were calculated using the real HF/MC6-311G** frequencies and, in addition. the HF/3-21G** mode at 81 cm-'. @

HF/3-21G** frequency and the imaginary HF/MC6-311G** frequency are very similar wagging modes of the bridging oxygen. In light of the many low frequencies, it is very difficult to state with conviction that any of these structures are or are not true minima. This uncertainty is reemphasized when electron correlation effects are included in the optimizations: Conformer C, the lowest-energy structure that we could locate on the HF potential energy surface, is no longer a stationary point for MP2/MC6311G** optimization. The effect of correlation is to strengthen the weaker H bridge so that conformer C collapses to the doublyH-bridged structure (conformer B). Conformer B remains 2.4 kcal mol-' lower in MP2 energy than the singly-bridged structure A (Table 5). With our computer hardware, we were unable to calculate MP2 frequencies for disilicic acid in order to verify the MP2 structures as minima. Thus, in order to estimate free-energy differences between the conformers, we used the real HF/MC6-311G** frequencies and, in addition, the HF/3-21G** mode at 81 cm-' for conformer B. Upon inclusion of zero-point energy, enthalpic, and entropic corrections, the free energy of conformer B is only 0.4 kcal mol-' below that of conformer A. For all three conformers, the use of MC6-311G** caused bond lengths to shrink compared to previous results,6,17,99,1M) as typically occurs with HF calculations when the quality of the basis set improves. As the Si-0 bonds shrink, the Si0-Si angle increases, as discussed above for HsSiOSiH3. Thus, the HF/MC6-311G** geometries have larger Si-0-Si angles than previous results for these molecular models. The MP2 corrections increased bond lengths by 1-3 pm compared to their values in the HF/MC6-311G** structures (Table 6). Note that several of the larger changes due to inclusion of electron correlation involved H-bonded atoms. The donor 09-H13 bond length in conformer A increased by 2.1 pm, while the receptor Si'-07 bond length increased by 2.5 pm: Meanwhile, the H-bridge length decreased by a full 62 pm, from 2.77 to 2.15 A! The increased strength of the

Silicate Molecular Models

J. Phys. Chem., Vol. 98, No. 48, I994 12551

TABLE 6: Optimized Geometric of H&iZO,, Using the MC6-311G** Basis Set

conformer A propertyb

HF

MP2

1.609 1.618 1.618 1.623 1.627 1.618 1.624 1.618 0.939 0.939 0.940 0.940 0.939 0.94 1 145.8 113.8 109.7 104.9 106.3 109.5 111.8 1 19.4 120.0 119.8 119.5 120.0 119.3 148.0 -89.9 29.8 116.3 - 124.4 -1.6 89.5 33.5 165.9 -166.1 -46.9 -52.7

1.629 1.647 1.636 1.640 1.652 1.638 1.645 1.633 0.957 0.957 0.958 0.957 0.957 0.962 134.4 115.4 109.2 104.1 104.8 109.7 111.6 115.0 114.7 115.8 114.7 114.9 112.9 153.7 -84.7 35.0 105.6 -134.2 -11.1 86.8 31.0 163.9 -166.3 -52.1 -39.7

conformer HF MP2 1.615 1.615 1.628 1.614 1.623 1.614 1.623 1.627 0.939 0.940 0.939 0.942 0.942 0.939 139.3 104.2 113.5 109.6 113.4 109.5 104.2 120.2 119.7 118.3 119.7 118.2 120.2 169.3 -70.3 51.3 50.1 168.1 -71.5 172.2 83.8 32.9 83.7 32.9 172.2

1.640 1.640 1.653 1.631 1.639 1.631 1.639 1.653 0.957 0.957 0.957 0.963 0.963 0.957 125.6 102.5 113.8 109.8 113.7 109.8 102.5 116.5 115.0 112.0 115.0 112.0 116.5 169.8 -68.6 52.8 52.4 169.4 -69.0 178.7 82.7 33.1 82.8 33.1 178.8

0.004

conformer c

HF

1.609 1.618 1.630 1.616 1.622 1.617 1.623 1.621 0.939 0.940 0.939 0.940 0.943 0.939 144.5 104.2 114.0 109.9 113.0 108.8 105.5 120.1 119.8 1 19.5 119.3 118.0 119.7 - 154.4 -34.9 87.5 13.0 130.4 -109.2 170.5 85.4 32.1 84.2 31.0 168.8

Bond lengths in angstroms; bond angles and torsions in degrees. Atom numberings taken from Figure 1.

hydrogen bond due to electron correlation contributed to an 11.6" decrease in the Si-0-Si angle when correlation was "turned on". Similarly, the two H-bridge lengths of conformer B decreased by 50 and 51 pm when electron correlation effects were accounted for. The resultant Si-0-Si angle of conformer B was 13.7" smaller at the MP2MC6-3 11G** level as compared to its HF-optimized value using the same basis set. For each conformer, the Si-0-H bond angles decrease by 4-6" when their determination includes electron correlation corrections. One would expect correlation to influence torsional angles e ~ p e c i a l l y , ~since , ~ ~ torsions are typically determined by nonbonded interactions. Indeed, for conformer A, the framework (Si-0-Si-0) torsions changed by up to 10.7". Averaging the framework torsions in Table 6 indicates that inclusion of electron correlation caused the Si, tetrahedron to rotate by a net 4.6" relative to the Si3 tetrahedron. Torsional angles were less affected by MP2 corrections for conformer B; the two H bonds confer extra rigidity to the structure. In Figure 3, the total MP2 electron density for H6Si207 and its difference from the I-IF density are plotted. If the minimum in the electron density along the Si-0 bond is used as a measure of atomic size,101J02 the silicon "radius" is found to be 0.36 A, while that of oxygen is 1.28 A. These are reasonably close to, but several picometers below, the typical ionic radiiIo3-Io5 reported for Si and 0 in silicates. The bonding is seen to be

0.014

0.00 0.00 -0.002 -0.005

0.003

/ I

-.. .

,

.a003

Figure 3. Total MpuMC6-311G**electron density for &Si207 (top) and density difference map in the Si-0-Si plane (bottom) at the MP2optimized geometry for H6Si207. The difference map illustrates the effects of electron correlation on the electron density by subtracting the total HF density from the total MP2 density. All numerical density values are given in atomic units.

covalent even in the Si-0-Si plane. As discussed above in conjunction with Figure 2, Figure 3 shows that MP2 corrections enhance the covalency of Si-0 bonding by increasing the electron density in the lowest-density region along the Si-0 bond and in regions of diffuse density overlap. While inclusion of correlation would be expected to transfer electron density toward diffuse regions of the molecule, it may be significant that electron correlation improves orbital overlap at the weakest point in the Si-0 bond.

Discussion: Comparison of the Molecules as Models for Silicate Minerals The calculations reported above used four uniform levels of quantum mechanical rigor to examine the molecular and electronic structures of disiloxane, orthosilicic acid, and disilicic acid. Thus, we are afforded an opportunity to directly compare the ab initio predictions and perhaps to assess the adequacy of these three molecules to represent structural units of mineral silicates. Disiloxane and its monomeric analogue, silanol (SiH3OH), have a rich history 12,13.15,16,37,43,670,72,74,75,78,79,92.l~~ 123 as molecular models for silicate functional groups because they are small molecules. Indeed, MC6-311G** calculations on

12552 J. Phys. Chem., Vol. 98, No. 48, 1994

Teppen et al.

TABLE 7: Average Calculated Si-0 Bonded Distances in Three Molecules9 method HF/MC6-311G** MP2/MC6-311G** HF/MC6-311G(2d,2p) MP2/MC6-311G(2d,2p)

Si-Obr Bonds

Si-OH bonds

disiloxane &Si207

Si(0H)d )I6Si207

1.621 1.640 1.623 1.645

1.614 1.638 1.609 1.637b

1.622 1.641 1.614 1.634

1.622 1.641 1.614 1.637*

Oxygens (G) bridging two silicans are distinguished from hydroxyl oxygens (OH). Frozen-core optimization, of conformer B only, due to hardware limitations.

H3SiOSiH3, Si(OH)4, and H6Si207 use 169, 201, and 355 primitive Gaussian basis functions to represent 106, 122, and 2 14 molecular orbitals, respectively. Each single-point MP2 calculation of wave function and force for &Si207 required about 42 times the CPU time as the same calculation for H3SiOSiH3. Besides being computationally expensive, the former molecule is more difficult to optimize, so disiloxane and its analogues remain very currently attractive12313915, 16,72,74,92J18,119J21-I23 as models for the bridging oxygen in silicates, including zeolites. It is of scientific interest whether this attractiveness is jusifiable on grounds other than computational convenience. A postulate at the outset is that, of the three molecules studied in this paper, H6Si207 should most closely model the properties of mineral silicates. In H6Si207, both s i atoms are fully coordinated by 0, as they would be in minerals, and the bridging 0 is surrounded by two proper coordination shells. Calculations of the quality reported here have not been done for larger silicate fragments, so no comparisons can be drawn with computational results that would be expected to better describe silicate minerals. Comparisons, if any, of the H6Si207 ab initio results will have to be made with experimental data, despite the lack of precise correspondence between ground-state ab initio results and experimental data that are thermal averages. Table 7 presents Si-0 bond lengths for the three molecules: All predicted bond lengths at all levels of theory fall within the large range of experimentally determined Si-0 bond lengths found in silicates, 1.55- 1.70 8,.83,95-97,124 The average experimental Si-0 distance in silicates is 1.626 A,124 Note that, for Si(OH)4, the highest quality calculations seem to be converging toward that value. That is, either an MP2 calculation with a basis set larger than MC6-311G(2d,2p) or a better treatment of correlation effects using the MC6-311G(2d,2p) basis should give average Si-0 bond lengths below 1.634 A, especially since there is no indication that we have reached the basis set limit. Although &Si207 was too large for full MpuMC6-311G(2d,2p) calculations to be executed by our computers, the frozen-core MP2 optimization of conformer B causes the Si-0 bond lengths to decrease as well, but by less than the full MP2 calculation on Si(OH)4. Otherwise, at each level of theory, the average Si-0 bond lengths found for Si(OH)d and &Si207 agree exactly. For H3SiOSiH3, in contrast, the Si-0 bond length seems to be increasing as the quality of the calculation increases. This behavior is contrary to expectation and to the trends found in the recent systematic Hartree-Fock study of basis set effects on disiloxane geometries.12 However, as discussed above in the disiloxane section, the highest quality ab initio results would be expected to converge to an Si-0 bond length somewhat above the experimental valuea of 1.634 A. This would still be in agreement with the mean bond length found in mineral silicates. Other bonded interactions cannot be compared among all three molecules, since H3SiOSiH3 has no hydroxyl groups and

Si(OH)4 has no bridging 0. It is difficult to compare values for the Si-0-Si angle in H3SiOSiH3 and H&&, since the angle is so flexible and so strongly influenced by the hydrogen bonding in H6Si207. Nevertheless, as remarked earlier, trends of Si-0 bond length as a function of Si-0-Si angle have been found to be descriptive of silicate minerals.83 The dependence of the bond length on the angle in &Si207 has been calculated using the STO-3Gg3and 6-31G*36 basis sets. In the former case, the functional dependence was seen to be too steep compared to experimental results, while in the latter case agreement was better but the slope was too shallow. In neither of these cases were full optimizations done for all degrees of freedom besides the angle. Even so, Lasaga and Gibbs6 showed that the 6-31G* results reproduced the experimentallZ5compressibility of the Si-0-Si angle in quartz resonably well. Purton et al.126used density functional methods on the central nine atoms of a larger cluster in which all Si atoms were fully coordinated by 0. The potential energy functions they derived were able to very closely match other experimental data93for Si-0-Si as a function of pressure. Recent high-quality Hartree-Fock calculations12 of H3SiOSiH3 showed that the Si-0 bond length seems to increase much too slowly as Si-0-Si decreases, compared with mineral silicate data. The authors12 conjecture that the Si-0 force constant in disiloxane is too large to accurately model silicates because Si is not fully coordinated by 0. Their hypothesis is corroborated by the evid e n ~ e ~ , ~ , that ~ - Si-0-Si ~ ~ , ~ ~force , ~ fields ~ . ~ derived ~ ~ from computations on H&i2O7 and larger fragments are very similar to those found in minerals. A comparison of hydroxyl parameters in Tables 4 and 6 indicates that, at the same level of theory, the 0-H bond lengths and the Si-0-H angles are very similar for both Si(OH)4 and H6Si~07. Such a close correspondence between bond lengths and angles for both Si-0 and 0-H interactions in these two molecules again indicates that predicted geometries of structures containing Si tetrahedrally coordinated by oxygen are likely to scale properly as tetrahedra are linked. Of course, proper calculation of nonbonded interactions is also crucial when the ultimate goal is to model dynamics of the silicate lattice or interactions of the lattice with adsorbates. Nonbonded Si-Si distances predicted for H3SiOSiH3 and HgSi207 are compared in Table 8. In every case, inclusion of correlation decreases the separation of the two Si. Experimental Si-Si separations for silicate minerals range from 2.9 to 3.3 A,83with that of a-quartz being nearly average at 3.058 A.93,94 Again, this is necessarily a vibrationally averaged distance that decreases slightly to 3.053 8, at 13 K.94 The ab initio distances would thus be expected to be slightly shorter than the latter value. The Si-Si separations for both H3SiOSiH3 and H6Si207 are in decent agreement with the experimental data, although those for H&07 tend to be shorter, probably due to intertetrahedral H bonds. The ab initio result for H3SiOSiH3 agrees rather well with the experimentala Si-Si distance of 3.107 8, for disiloxane. The Si(OH)4 and H6Si207 results may be further compared by examining the nonbonded distances between geminal (13) oxygens, as listed in Table 8. Again, the agreement is excellent (within 0.003 A at worst) when the geometries are compared at the same level of theory. The 0-0distance seems to be converging toward about 2.66 A, which is slightly higher than the experimental average of 2.63 A for quartz.93 Note that inclusion of correlation lengthens the 1,3-0-0 nonbonded distance, as a result of lengthening all Si-0 bonds. Based on these changes in nonbonded distances, the effects of correlation on a model crystal structure can be imagined.

J. Phys. Chem., Vol. 98, No. 48, 1994 12553

Silicate Molecular Models TABLE 8: Selected Nonbonded Distance Calculated for the Silicate Molecules nonbonded interaction

molecule

I’

basis set MC6-31 1G(2d.2~) MC6-311G** MC6-311G** MC6-31lG(2d.2~) MC6-311G** MC6-311G** MC6-311G** MC6-31lG(2d,2p) MC6-311G** MC6-311G** MC6-31lG(2d.2~) MC6-311G** MC6-31 1G** MC6-31lG(2d.2~) MC6-311G** MC6-311G** MC6-311G** MC6-311G(2d92p) MC6-311G**

nonbonded dist, A HF MP2 3.111 3.085 3.085 3.020 3.028 2.918 2.904 2.963 3.073 2.153 2.769 2.781 2.269 2.249 2.632 2.444 2.679 2.644 2.635 2.668 2.677 2.643 2.678 2.644 2.672 2.632 2.644 2.916 3.290 3.424 3.045 3.295 3.031 3.153

-0.008

,

-0.018

\,

/‘

\

4.028

‘ \,

’

0.003

“When more than one instance of the given interaction occurs in the molecule, the distances are averaged. The 1,3-0-0 distance is thus the average of six such distances in each tetrahedron. The 0-H and 1,5-0-0 distances are for only those atoms in the hydroxyl groups that are involved in the hydrogen bridges of H6Si207.

MP2-optimized structures result in larger individual tetrahedra relative to HF-optimized structures, since the 1,3-0-0 distance increases. However, individual tetrahedra will pack closer together due to the dispersion forces that are included in the correlated calculations, as evidenced by the decrease in Si-Si. Thus, a “correlated” quartz crystal would be more dense than a “Hartree-Fock” quartz crystal, even though each tetrahedron in the former would be larger. The fact that Si-0-Si optimizes to reasonable silicate values in HaSi20-1at the HF/MC6-311G** level of theory, whereas the same calculation on disiloxane produces a linear structure, indicates that the electronic environment of oxygen may be quite different in the two situations. However, comparison of MP2 electron densities near the bridging oxygen for H3SiOSiH3 and H&07 (Figures 2 and 3) shows no obvious differences. Indeed, in comparing the electron density at a point 1.O A from 0, the density is 0.041 au for HsSiOSiH3 and it is only slightly higher at 0.044 au for H6Si207. A major difference between the two cases is the electrostatic potential that results from integrating the MP2 electron densities, as displayed in Figure 4. For disiloxane, the minimum in-plane molecular electrostatic potential is -0.0387 am., with a mean potential gradient of 0.0098 au. In contrast, the molecular electrostatic potential of H&07 exhibits a potential minimum of -0.0521 au and a gradient Qf -0.0179 au leading into that minimum. In addition, the potential minimum is about 0.1 A closer to the oxygen of HaSi20-1than that of H3SiOSiH3. Thus, a positive charge would feel a stronger force of attraction as well as find a lower-energy equilibrium near &Si207 than it would near H3SiOSiH3. The electrostatic potential plots (Figure 4) are consistent with our calculations of ab initio partial charges on the bridging Obr in H3SiOSiH3 and H&07 (Table 9) that show Obr in H&07 to have an effective negative charge about 0.24 electrons greater in magnitude than the Ob in H3SiOSiH3. This is a huge difference, since electrostatics would be expected to contribute most of the interaction energy between Q r and polar or ionic species. This charge difference is further evidence that &Si207 should be a much more successful molecular model for silicate minerals than H3SiOSiH3.

Figure 4. Molecular electrostatic potential maps for H3SiOSiH3 (top) and H6Si207 (bottom). The potential is given in atomic units.

Table 9 illustrates the importance of the method chosen for estimating atomic partial charges. The standard Mulliken has long been recognized as a crude estimate of the total charge on an atom, and it can be seen from Table 9 that the Mulliken charges on all atoms except hydrogen increase dramatically as the size of the basis set increases. The CHELP method of Chirlian and Francl131 is much more physically meaningful, assigning charges by integrating the molecular electrostatic potential at a grid of points and choosing atomic charges that best reproduce the grid. The CHELPG method132 improved upon CHELP by using a finer and more regularly spaced grid and is expected to yield the most meaningful results of the three methods. Note that, if the Mulliken charges were used to model the electrostatic potential near Obr, then the potential energy well would be deeper near H3SiOSiH3 than near H6Si207, which is opposite to the results from integrating the electron density (Figure 4). The CHELPG results exhibit very consistent trends across basis set and level of quantum theory. For every atom, either increasing the size of the basis set or adding electron correlation effects decreased the magnitude of the atom’s partial charge (with the exception of H3SiOSiH3 with the MC6-311G** basis set, which predicted a linear molecule at the HF level but a bent molecule at the MP2 level). Adding a second set of polarization functions increases charge density in the bonding regions of a molecule,81 thus lowering the absolute value of

12554 J. Phys. Chem., Vol. 98, No. 48, 1994

Teppen et al.

TABLE 9: Effect of the Correlation on ab Znifio Partial Charges: Calculated by Two Different Methods, for the Atoms in Three Silicate Molecules

basis set and charge calculation procedure MC6-311G** basis set MC6-31lG(2d.2~)basis set type of calcn Multiken charge CHELPG charge Mulliien charge CHELPG charge atomband molecule Si in HlSiOSiH? HF 1.151 0.890 1.434 0.860 MP2 1.150 0.912 1.205 0.815 Si in &Si207 HF 1.676 1.786 2.475 1.551 MP2 1.352 1.515 HF Si in Si(OH)4 1.611 2.412 1.546 1.763 ME 1.306 1.507 1.943 1.606 G in HrSiOSiHs HF -0.977 -0.565 -1.323 -0.550 MP2 -0.988 -0.595 -1.104 -0.528 -0.941 Obr in &Si207 HF -0.922 -1.332 -0.794 MP2 -0.790 -0.788 0” in H6Si207 HF -0.712 -0.888 -0.854 -0.934 MP2 -0.613 -0.833 OHin Si(OH)4 HF -0.881 -0.858 -0.71 1 -0.935 ME -0.706 -0.840 -0.761 -0.771 HOin H6Si~07 HF 0.310 0.492 0.285 0.469 MP2 0.294 0.459 HF 0.308 HOin Si(0H)d 0.278 0.471 0.494 MP2 0.275 0.445 0.305 0.463 H in H3SiOSiH3 HF -0.221 -0.203 -0.258 -0.195 MP2 -0.219 -0.205 -0.218 -0.184 Charges in electron charges. When there are multiple atoms of a given type in a molecule, their charges are averaged. For &Si207 with the MC6-311G** basis set, charges on the atoms of the three conformers were also averaged together (two conformers for MP2 since conformer C was not stable when electron correlation was included in the calculations). &Si207 could not be optimized at the full MP2/MC6-311G(2d,2p)level due to computer memory limitations, and the MP2 frozen core partial charges for &Si207 are not comparable to the full MP2 partial charges obtained for the other systems. G refers to oxygen bridging two silicons, OH refers to hydroxyl oxygen, and Ho refers to hydroxyl hydrogen. each effective atomic charge. Electron correlation, on the other hand, removes charge density from the bonding regions and redistributes it around the edges of a molecule, again reducing the magnitude of each effective point charge.62 The changes in electron density near the Si-0-Si functional group that lead to this reduction in the effective charge on %when correlation effects are accounted for are displayed in Figures 2 and 3. The pattem is very similar for the two rather dissimilar molecules H3SiOSiH3 and H6Si207. In addition to the CHELFG charge on & in disiloxane being far less negative than in H6Si207, the charge on Si in H3SiOSiH3 is far less positive, indicating that both 0 and Si in H3SiOSiH3 would be expected to interact too weakly with the environment compared to 0 and Si in silicate minerals. Sauer and Ahlrichs116calculated the partial charges for the atoms of silanol including both a large basis set and correlation effects: The magnitudes of those charges were even less than those of disiloxane, with Si and 0 possessing 0.42 and -0.35 charge, respectively. Further examination of Table 9 shows that the CHELPG charges on corresponding atoms in Si(OH)4 and H6Si207 at a given level of theory are almost exactly the same in every instance. Thus, electrostatic potentials, partial atomic charges, and, hence, reactivities of atoms in silicate fragments seem to be a function of the proper coordination of Si by four 0 atoms rather than a function of molecular fragment size. This would imply that Si(OH)4 would make an adequate model for interactions involving the silicate hydroxyl but that silanol would not and that H6Si207 would provide the smallest quantitative model for the bridging 0 in silicate minerals. The differences among these molecules can be used to illuminate the results of several studies that sought to examine the properties of solid silicates by using molecular models. The dissimilarities between silanol or disiloxane models, on the one hand, and silicic acid or disilicic acid models, on the other, result in major discrepancies between predicted chemical behavior of the two types of models.

For example, Ugliengo and c o - w o r k e r ~ ~13s1~ l4 , ~used ~ J the disiloxane analogue silanol to model the surface hydroxyl of silica. Even high-quality MP2 calculations of silanol interactions with H207*and NH379gave equilibrium constants that were 2 orders of magnitude too low for the adsorption of these two species onto silica. The a ~ t h o r sattributed ~ ~ , ~ ~these poor results to their model’s inability to account for entropic contributions from perturbed crystal vibrations. An alternative and simpler explanation, based on the present study, seems to be that silanol provides a poor electrostatic model for a silicate hydroxyl. Another recent study16 showed that ammoniation energies calculated using disiloxane and trisiloxane analogues are much too low compared to experimental values for zeolites. Molecular electrostatic potentials indicate that adsorption energies based on silicic acid or disilicic acid models for the silica hydroxyl would be much more negative and thus closer to the observed values. Brand et al.llg calculated proton affinities (PA’s) of a bridging o b r atom (atom 0(24)133 of zeolite ZSM-5) for a series of clusters of increasing size and found the PA’s to vary by a rather large 44 kcal mol-’, without seeming to converge. However, most of the variation was between PA’s for clusters containing disiloxane-type groups (H3SiObrAlH3- and (H3Si0)3SiObrA1(OSiH3)3-) and clusters in which Si and A1 were fully coordinated by four oxygens ((H0)3SiObrA1(OH)3- and ((H0)3SiO)3Si@,Al(OSi(OH)3)3-). The PA’s for the two HOterminated clusters were similar and generally highest and thus agreed best with experimental values for the ze01ite.l~~Furthermore, other calculations135using the same basis set and a ring of three Si and one A1 saturated by 0 yielded a PA just between those of the smaller and larger 0-saturated clusters of the previously discussed study.118 Brand and associates’l8 tentatively explained the trends through electrostatic effects: Adding a coordination sphere of positively charged Si decreases PA, while a sphere of negative 0 increases PA. The data displayed in Figure 4 provide quantitative support for their qualitative explanation. Our data

Silicate Molecular Models

J. Phys. Chem., Vol. 98, No. 48, 1994 12555

show that SiH3 rather than Si(OH)3 substituents on Obr result polarized basis sets should be used when possible. The 6-31G* in a 26% decline in the depth of the nearby electrostatic potential and Dunning (9sSp) basis sets seem too small, causing MP2 well. Brand’s118 data imply that -SiH3 groups even several methods to overcorrect rather badly for correlation. The MC6coordination spheres away from the central &r cause major 311G** basis performed adequately but was clearly short of changes in the electrostatic potential near Obr. The adsorption the basis set limit. Adding a second set of polarization functions of ammonium ion by the bridging 0 was also studied by Brand to form MC6-311G(2d,2p) decreased bond lengths by about et a1.118 The adsorption energy followed the same trends as 0.8 pm in both HF and MP2 calculations. Based on a those for PA, with the two HO-terminated fragments giving comparison of the present study with previous results for comparable results that showed greater affinity for than disiloxane,12the Dunning (12s,9p) basis augmented with (2d,either of the two HsSi-terminated clusters. In a subsequent 2p) polarization functions should give very similar geometries paper, Brand et al.lZ3argue that OH groups should not be used and slightly lower energies than MC6-311G(2d,2p) in silicate to terminate molecular models for silicates because the eleccalculations. trostatic potential minimum near the hydroxyl oxygen is lower (d) When partial atomic charges are desired for modeling than that near the bridging oxygen of H3SiOSiHzOH. However, purposes, the CHELPG132method gave far more consistent we have shown (Figure 4) that the OH termini also contribute results than the standard Mulliken population analysis and to a lower electrostatic potential near the bridging oxygen, so differentiated much more effectively between atoms in the Brand’s datalz3may also be used to arguefor the use of clusmolecular models considered here. ters terminated by OH groups. Indeed, their own earlier study118 The computational demands of the techniques outlined above indicated that bridging oxygens of HO-terminated moleare considerable. However, they seem necessary to the develcules exhibited proton affinities that agreed better with experiopment of accurate ab initio force fields since geometries, partial ment. charges, and energies are so dependent on correlation effects, Both sets of s t ~ d i e s ~ ~discussed , ~ ~ . ~above ~ , ~show ~ ~that , ~ ~ ~which in turn are sensitive to the quality of the basis set. oxygen atoms in molecular clusters terminated by SiH3 interact Acknowledgment. We are grateful for support by the far too weakly with polar molecules and ionic species, compared University of Arkansas Graduate School, the College of to the analogous 0 in silicate minerals. The bridging 0 in HOAgriculture and Home Economics, and USDNCSRS National terminated fragments interacts more strongly with these species.118,135 Research Initiative Grant 93-37102-9557. These results can be explained on the simple basis of Supplementary Material Available: One table listing the electrostatic potentials (Figure 4, Table 9) near the various optimized geometries from the ab initio study of H6Si207. molecular models. Since molecular electrostatic potentials and conformer B (1 page). Ordering information is given on any the partial charges fit to them are the most important quantum current masthead page. chemical and molecular mechanics factors that influence interactions with ionic or polar species, Si(OH)4 and (HO)3SiOSiReferences and Notes (OH)3 should be far superior to H3SiOH and H3SiOSiH3 as models for functional groups in silicate minerals. (1) Clementi, E. Computational aspects for large chemical systems;

m+

Conclusions The geometries of three molecules often used to represent structural units of silicate minerals were fully optimized to examine the effects of basis set size and inclusion of electron correlation. Based on these results, some general recommendations may be given for ab initio calculations of molecules as models for exploring the properties of silicate minerals and their interactions with adsorbates. (a) The molecular fragment chosen should at least contain Si atoms (and by extension other metals) that are properly coordinated by oxygen. This advice seems trite but needs to be said in view of the recent upsurge in paper~12,13,15.16,31,72,74,78,92,112,119,121-123~136 representing silicate structures by silanol and disiloxane analogues. Fortunately, geometries and charges calculated for Si(OH)4 agree very well with those for H6Si207, adding further evidence to arguments that local electronic effects are the most important determinant of crystal geometries. (b) The extra effort needed to include correlation effects would seem to be worthwhile, in view of the importance of accurate local electronic structure and geometry. The MP2 corrections employed in this work tended to lengthen all bonds and reduce all partial charges systematically,but had sometimes large and often less systematic effects on bond angles, torsions, and nonbonded distances. In addition, the global minimum of H6Si207 according to Ell3 calculations seems to not even be a stationary point when electron correlation is included in the optimization. (c) This research supports and extends to other silicate molecules previous results, 12,13 indicating that large, multiply

Lecture Notes in Chemistry 19; Springer-Verlag: Berlin, 1980. (2) O’Keeffe, M.; Newton, M. D.; Gibbs, G. V. Phys. Chem. Miner. 1980,6, 305-312. (3) O’Keeffe, M.; McMillan, P. F. J. Phys. Chem. 1986, 90, 541542. (4) Hess, A. C.; McMillan, P. F.; O’Keeffe, M. J. Phys. Chem. 1987, 91, 1395-1396. (5) Lasaga, A. C.; Gibbs, G. V. Phys. Chem. Miner. 1987, 14, 107117. (6) Lasaga, A. C.; Gibbs, G. V. Phys. Chem. Miner. 1991, 17, 485491. (7) Tsuneyuki, S.; Tsukada, M.; Aoki, H.; Matsui, Y. Phys. Rev. Lerr. 1988, 61, 869-872. (8) Catlow, C. R. A.; Price, G. D. Nature 1990, 347, 243-248. (9) McMillan, P. F.; Hess, A. C. Phys. Chem. Miner. 1990, 17, 97107. (10) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Phys. Rev. Lett. 1990, 16, 1955-1958. (11) Kramer, G . J.; Farragher, N. P.; van Beest, B. H. W.; van Santen, R. A. Phys. Rev. B 1991, 43, 5068-5080. (12) Nicholas, I. B.; Winans, R. E.; Harrison, R. J.; Iton, L. E.; Curtiss, L. A.; Hopfinger, A. J. J. Phys. Chem. 1992, 96, 7958-7965. (13) Nicholas, I. B.; Winans, R. E.; Harrison, R. J.; Iton, L. E.; Curtiss, L. A.; Hopfinger, A. J. J. Phys. Chem. 1992, 96, 10247-10257. (14) Purton, J.; Jones, R.; Heggie, M.; Oberg, S.; Catlow, C. R. A. Phys. Chem. Miner. 1992, 18, 389-392. (15) Teunissen, E. H.; van Santen, R. A,; Jansen, A. P. J.; van Duijneveldt, F. B. J. Phys. Chem. 1993, 97, 203-210. (16) Kassab, E.; Fouquet, J.; Allavena, M.; Evleth, E. M. J. Phys. Chem. 1993, 97, 9034-9039. (17) Hill, J.-R.; Sauer, J. J. Phys. Chem. 1994, 98, 1238-1244. (18) Pisani, C . ; Dovesi, R. Znr. J. Quantum Chem. 1980, 17, 501-516. (19) Pisani, C.; Dovesi, R.; Roetti, C. Hartree-Fock ab initio treatment of crystalline systems; Lecture Notes in Chemistry 48; Springer-Verlag: Berlin, 1988. (20) Dovesi, R.; Pisani, C.; Roetti, C.; Causa, M.; Saunders, V. R. CRYSTAL88. An ab initio all-electron LCAO-Harrree-Fock program for periodic solids. Program no. 577; QCPE: Bloomington, IN, 1988. (21) Dovesi, R.; Pisani, C.; Roetti, C.; Silvi, B. J. Chem. Phys. 1987, 86, 6967-6911.

12556 J. Phys. Chem., Vol. 98, No. 48, 1994 (22) Dovesi, R.; Roetti, C.; Freyria-Fava, C.; Aprh, E.; Saunders, V. R.; Harrison, N. M. PhiIos. Trans. R. SOC.London, Ser. A 1992,341,203210. (23) Hess, A. C.; Saunders, V. R. J. Phys. Chem. 1992, 96,43674374. (24) Hess, A. C.; McCarthy, M. I.; White, J. C.; Anchell, J.; Nicholas, J. B.; Thompson, M. R. Chem. Des. Autom. News 1993,8, 1,36-43. (25) White, J. C.; Hess, A. C. J. Phys. Chem. 1993,97,8703-8706. (26) White, J. C.; Hess, A. C. J. Phys. Chem. 1993,97,6398-6404. (27) Tossell, J. A.; Gibbs, G. V. Phys. Chem. Miner. 1977,2,21-57. (28) Gibbs, G. V. Am. Mineral. 1982,67,421-450. (29) Sauer, J.; Zahradnik, R. lnt. J. Quantum Chem. 1984,26,793822. (30) Sauer, J. Chem. Rev. 1989,89, 199-255. (31) Lasaga, A. C. Rev. Mineral. 1990,23, 17-85. (32) Tossell, J. A.; Vaughan, D. J. Theoretical geochemistry: Application of quantum mechanics in the earth and mineral sciences; Oxford University: New York, 1992. (33) Sauer, J. In Modelling of structure and reactivity in zeolites; Catlow, C. R. A., Ed.; Academic: London, 1992; pp 183-216. (34) Burdett, J. K.; Caneva, D. C. Inorg. Chem. 1985,24,3866-3873. (35) Newton, M. D.; Gibbs, G.V. Phys. Chem. Miner. 1980,6,221246. (36) O’Keeffe, M.; Domengks, B.; Gibbs, G. V. J. Phys. Chem. 1985, 89,2304-2309. (37) Gibbs, G. V.; D’Arco, P.; Boisen, M. B., Jr. J. Phys. Chem. 1987, 91,5347-5354. (38) Navrotsky, A.; Geisinger, K. L.; McMillan, P.; Gibbs, G.V. Phys. Chem. Miner. 1985,11, 284-298. (39) Kramer, G.J.; de Man, A. J. M.; van Santen, R. A. J. Am. Chem. SOC. 1991,113,6435-6441. (40) Frey, R. F.; Coffin, J.; Newton, S. Q.;Ramek, M.; Cheng, V. K. W.; Momany, F. A.; Schger, L. J. Am. Chem. SOC. 1992, 114, 53695377. (41) Momany, F. A.; Rone, R.; Frey, R. F.; Schger, L. Chem. Des. Autom. News 1992,7, 1. (42) Momany, F. A.; Rone, R.; Frey, R. F.; S c h ~ e r L. , QUANTA Parameter Handbook, Release 3.3;Molecular Simulations: Burlington, MA, 1992. (43) Sauer, J. J. Phys. Chem. 1987,91,2315-2319. (44) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A,; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian 92,Revision B; Gaussian: Pittsburgh, PA, 1992. (45) Sellers. H. L.: Klimkowski. V. J.: Schtifer. L. Chem. Phvs. Lett. 1978,58,541-544. ’ (46) Pulay, P. Mol. Phys. 1969,17, 197-204. (47) Pulay, P.; Fogarasi, G.; Pang, F.; Boggs, J. E. J. Am. Chem. SOC. 1979,101,2550-2560. (48) Chsky, P.; Hess, B. A., Jr.; Schaad, L. J. J. Am. Chem. SOC. 1983, 105,396-402. (49) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72,650-654. (50) McLean, A. D.; Chandler, G.S. J. Chem. Phys. 1980,72,56395648. (51) Wong, M. W.; Gill, P. M. W.; Nobes, R. H.; Radom, L. J. Phys. Chem. 1988,92,4875-4880. (52) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 1982,77,3654-3665. (53) Raffenetti, R. C. J. Chem. Phys. 1973,58, 4452-4458. (54) Binkley, J. S.; Pople, J. A. lnr. J. Quantum Chem. 1975,9,229236. (55) Frisch, M. J.; Pople, J. A.; Binkley, J. S. J. Chem. Phys. 1984,80, 3265-3269. (56) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem Phys. 1972,56, 2257-2261. (57) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973,28, 213222. (58) M~ller,C.; Plesset, M. S. Phys. Rev. 1934,46,618-622. (59) Pople, I. A.; Binkley, J. S.; Seeger, R. lnt. J. Quantum Chem. Symp. 1976,IO, 1-19. (60) Levine, I. N. Quantum chemistry, 4th ed.; Prentice Hall: Englewood Cliffs, NJ, 1991. (61) Wolinski, K.; Pulay, P. J. Chem. Phys. 1989,90, 3647-3659. (62) Wiberg, K. B.; Hadad, C. M.; LePage, T. J.; Breneman, C. M.; Frisch, M. J. J. Phys. Chem. 1992,96,671-679. (63) Handy, N. C.; Schaefer, H. F., III.J. Chem. Phys. 1984,81,50315033. (64) Almenningen, A.; Bastiansen, 0.; Ewing, V.; Hedberg, K.; Traetteberg, M. Acta Chem. Scand. 1963,17,2455-2460. (65) Barrow, M. J.; Ebsworth, E. A. B.; Harding, M. M. Acta Crystallogr. 1979,B35, 2093-2099. (66) Sauer, J.; Zurawski, B. Chem. Phys. Lett. 1979,65,587-591.

Teppen et al. (67) Grigoras, S.; Lane, T. H. J. Comput. Chem. 1987,8,84-93. (68) Abraham, R. J.; Grant, G.H. J. Computer-Aided Mol. Des. 1988, 2, 267-280. (69) Koput, J. Chem. Phys. 1990,148,299-308. (70) Earley, C. W. J. Comput. Chem. 1993,14, 216-225. (71) Dunning, T. H., Jr. J. Chem. Phys. 1971,55, 716-723. (72) Luke, B. T. J. Phys. Chem. 1993,97,7505-7510. (73) Ahlrichs, R.; S c h d , P.; Ehrhardt, C. J. Chem. Phys. 1985,82,890898. (74) Curtiss, L. A.; Brand, H.; Nicholas, J. B.; Iton, L. E. Chem. Phys. Lett. 1991, 184,215-220. (75) Lindsay, C. G.; Tossell, J. A. Phys. Chem. Miner. 1991,18,191198. (76) Shambayati, S.; Blake, J. F.; Wierschke, S. G.;Jorgensen, W. L.; Schreiber, S. L. J. Am. Chem. SOC.1990,112,697-703. (77) Dung, J. R.; Flanagan, M. J.; Kalasinsky, V. F. J. Chem. Phys. 1977,66,2775-2785. (78) Ugliengo, P.; Saunders, V.; Garrone, E. J. Phys. Chem. 1990,94, 2260-2267. (79) Ugliengo, P.; Saunders, V. R.; Gamone, E. Su@ace Sci. 1989,224, 498-514. (80) Magnusson, E. J. Comput. Chem. 1993,14, 54-66. (81) Magnusson, E. J. Am. Chem. SOC. 1993,115,1051-1061. (82) Gibbs, G. V.; Hamil, M. M.; Louisnathan, S. J.; Bartell, L. S.; Yow, H. Am. Mineral. 1972,57, 1578-1613. (83) Hill, R. J.; Gibbs, G. V. Acta Crystallogr. 1979,B35, 25-30. (84) Gibbs, G. V.; Meagher, E. P.; Newton, M. D.; Swanson, D. K. In Structure and bonding in crystals; O’Keeffe, M., Navrotsky, A., Eds.; Academic: New York, 1981; Vol. I, Chapter 9, pp 195-225. (85) Sauer, J. Chem. Phys. Lett. 1983,97,275-278. (86) De Almeida, W. B.; O’Malley, P. J. Chem. Phys. Lett. 1991,178, 483-487. (87) Ferris, K. F. J. Mol. Struct. 1992,257,499-506. (88) Hess, A. C.; McMillan, P. F.; O’Keeffe, M. J. Phys. Chem. 1988, 92,1785-1791. (89) Magnusson, E. J. Am. Chem. SOC. 1990,112,7940-7951. (90) De Almeida, W. B.; O’Malley, P. J. THEOCHEM 1992,257,305312. (91) De Almeida, W. B.; O’Malley, P. J. THEOCHEM 1992,258,379387. (92) De Almeida, W. B.; O’Malley, P. J. J. Chem. SOC., Faraday Trans. 1993,89,983-989. (93) Levien, L.; Prewitt, C. T.; Weidner, D. J. Am. Mineral. 1980,65, 920-930. (94) Lager, G. A.; Jorgensen, J. D.; Rotella, F. J. J. Appl. Phys. 1982, 53,6751-6756. (95) Ribbe, P. H., Ed. Orthosilicates, 2nd ed.; Rev. Mineral. 5; Mineralogical Society of America: Washington, DC, 1982. (96) Bailey, S. W., Ed. Hydrous phyllosilicates (exclusive of micas); Rev. Mineral. 19; Mineralogical Society of America: Washington, DC, 1988. (97) Griffen, D. T. Silicate crystal chemistry; Oxford University: New York, 1992. (98) Lasaga, A. C.; Gibbs, G.V. Phys. Chem. Miner. 1988,16,2941. (99) Burkhard, D. J. M.; De Jong, B. H. W. S.; Meyer, A. J. H. M.; van Lenthe, J. H. Geochim. Cosmochim. Acta 1991,55, 3453-3458. (100) Kubicki, J. D.; Sykes, D. Am. Mineral. 1993, 78,253-259. (101) Gibbs, G. V.; Spackman, M. A,; Boisen, M. B., Jr. A m Mineral. 1992,77,741-750. (102) Feth, S.; Gibbs, G. V.; Boisen, M. B., Jr.; Myers, R. H. J. Phys. Chem. 1993,97, 11445-11450. (103) Goldschmidt, V. M. Natunvissenschaffen 1926,14,477-485. (104) Pauling, L. J. Am. Chem. SOC. 1927,49,765-790. (105) Klein, C.; Hurlbut, C. S., Jr. Manual of mineralogy ( a f e r James D.Dana), 20th ed.; John Wiley and Sons: New York, 1985. (106) Sauer, J.; Hobza, P.; Zahradnik, R. J. Phys. Chem. 1980, 84, 3318-3326. (107) Hobza, P.; Sauer, J.; Morgeneyer, C.; Hurych, J.; Zahradnik, R. J. Phys. Chem. 1981,85, 4061-4067. (108) Mortier, W. J.; Sauer, J.; Lercher, J. A.; Noller, H. J. Phys. Chem. 1984,88,905-912. (109) Geerlings, P.; Tariel, N.; Botrel, A.; Lissillour, R.; Mortier, W. J. J. Phys. Chem. 1984,88, 5752-5759. (110) Chakoumakos,B. C.; Gibbs, G.V. J. Phys. Chem. 1986,90,996998. (111) Beran, S. J. Phys. Chem. 1990,94,335-337. (112) Lasaga, A. C.; Gibbs, G. V. Am. J. Sci. 1990,290,263-295. (113) Ugliengo, P.; Gamone, E. J. Mol. Catal. 1989,54,439-443. (114) Ugliengo, P.; Saunders, V.R.; Gamone, E. J. Phys. Chem. 1989, 93,5210-5215. (115) Gamone, E.; Ugliengo, P. Lnngmuir 1991, 7,1409-1412. (116) Sauer, J.; Ahlrichs, R. J. Chem. Phys. 1990,93,2575-2583. (117) Hill, J.-R.; Sauer, J.; Ahlrichs, R. Mol. Phys. 1991, 73, 335348.

Silicate Molecular Models (118) Brand, H. V.; Curtiss, L. A.; Iton, L. A. J. Phys. Chem. 1992,96, 7725-7732. (119) Bates, S.; Dwyer, J. J. Phys. Chem. 1993, 97, 5897-5900. (120) Langenaeker, W.; De Decker, M.; Geerlings, P.; Raeymaekers, P. THEOCHEM 1990, 207, 115-130. (121) Kramer, G. J.; van Santen, R. A.; Emeis, C. A.; Nowak, A. K. Nature 1993, 363, 529-531. (122) Jentys, A.; Grimes, R. W.; Gale, J. D.; Catlow, C. R. A. J. Phys. Chem. 1993,97, 13535-13538. (123) Brand, H. V.; Curtiss, L. A,; Iton, L. E. J. Phys. Chem. 1993, 97, 12773-12782. (124) Tsirelson, V. G.; Evdokimova, 0. A.; Belokoneva, E. L.; Umsov, V. S. Phys. Chem. Miner. 1990, 17, 275-292. (125) Hazen, R. M.; Finger, L. W.; Hemley, R. J.; Mao, H. K. Solid State Commun. 1989, 72, 507-511. (126) Purton, J.; Jones, R.; Catlow, C. R. A.; Leslie, M. Phys. Chem. Miner. 1993, 19, 392-400.

J. Phys. Chem., Vol. 98, No. 48, 1994 12557 (127) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833-1840. (128) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1841-1846. (129) Mulliken, R. S. J. Chem. Phys. 1955, 23, 2338-2342. (130) Mulliken, R. S. J. Chem. Phys. 1955, 23, 2343-2346. (131) Chirlian, L. E.; Francl, M. M. J. Comput. Chem. 1987,8,894-905. (132) Breneman, C. M.; Wiberg, K. B. J. Compuf. Chem. 1990, 11, 361-373. (133) Olson, D. H.; Kokotailo, G. T.; Lawton, S. L.; Meier, W. M. J. Phys. Chem. 1981, 85, 2238-2243. (134) Datka, J.; Boczan, M.; Rymarowicz, P. J. Catal. 1988,114,368376. (135) Kramer, G. J.; van Santen, R. A. J. Am. Chem. Soc. 1993, 115, 2887-2897. (136) Lasaga, A. C.; Gibbs, G. V. In Aquatic chemical kinetics: Reaction rates of processes in natural waters; Stumm, W., Ed.; John Wiley and Sons: New York, 1990; Chapter 9, pp 259-289.