10534
J. Phys. Chem. B 2000, 104, 10534-10542
Model Structures and Electron Density Distributions for the Silica Polymorph Coesite at Pressure: An Assessment of OO Bonded Interactions G. V. Gibbs Departments of Geological Sciences, Materials Science and Engineering and Mathematics, Virginia Tech, Blacksburg, Virginia 24061
M. B. Boisen, Jr.* Department of Mathematics, Virginia Tech, Blacksburg, Virginia 24061
K. M. Rosso Pacific NW National Laboratories, EnVironmental Molecular Sciences Laboratory, P.O. Box, K8-96, Richland, Washington 99352
D. M. Teter Geochemistry Department, Sandia National Laboratories, Albuquerque, New Mexico 87185-0750
M. S. T. Bukowinski Department of Geology and Geophysics, UniVersity of California at Berkeley, Berkeley, California 74720 ReceiVed: June 12, 2000; In Final Form: September 8, 2000
The crystal structure and the bond critical point, bcp, properties of the electron density distribution for the high-pressure silica polymorph coesite were generated for pressures up to ∼17 GPa, using first-principles calculations. The nonequivalent SiO bond lengths and the SiOSi and OSiO angles of the generated structures agree with those observed to within ∼1%. With compression, the SiO bond lengths and the variable SiOSi angles of the structures both decrease while the value of the electron density, F(rc), the curvatures, and the Laplacian of the electron density distribution at the bond critical points each increases slightly. As found in a recent modeling of the structure of low quartz, the calculated electron density distributions are nearly static and change relatively little with compression. The bcp properties of the model structure agree with those observed at ambient conditions to within ∼10%, on average, with several of the properties observed to correlate with the observed SiO bond lengths, R(SiO). This agreement is comparable with that observed for several other silicates. As predicted, the bonded radius of the oxide anion, the curvatures of F(rc) paralleling the bond paths and the Laplacian of F(rc) each correlates with the observed bond lengths. However, the observed F(rc) values and the curvatures of F(rc) perpendicular to the paths fail to show a correlation with the observed bond lengths. The ellipticity of the SiO bonds in both the model and the observed structures tends to decrease in value as the SiOSi angle approaches 180°, indicating that the bonds become more circular in cross sections as the angle widens. Ridges of electron density and bond critical points were found between the intertetrahedral oxide anions at each pressure. The existence of these features appears to be closely related to purely geometrical factors of the coesite structure rather than to bonded interactions. None of these features was found between the intratetrahedral oxide anions.
Introduction The electron density of a molecule or a crystal in a stationary state adopts a distribution wherein the total energy of the resulting configuration is minimized.1 A mapping of the distribution in a plane through the nuclei of a pair of bonded atoms generally displays two well-defined maxima connected by a relatively low-lying ridge of electron density. The two maxima define the positions of the atoms and the ridge follows the line in the distribution (the bond path) along which F(r) is a maximum relative to any neighboring line. The value of F(r) * To whom correspondence should be addressed. E-mail: boisen@ math.vt.edu.
typically decreases exponentially between the two maxima along the ridge to a point between the atoms where F(r) adopts a minimum value, and the gradient of the vector field, ∇F(r), of the electron density distribution vanishes. As ∇F(r) ) 0 at this point, the point is referred to as a bond critical point denoted rc. The distance between rc and the nuclei of each bonded atom is logically defined to be the bonded radius of the atom. The value of F(r) evaluated at rc has been found to serve as a measure of the strength of the interaction between a pair of bonded atoms, the greater the value of F(rc), the greater the strength and the shorter the bond. The bond critical point properties1 of a pair of bonded atoms including F(rc), the bonded radii together with the curvatures, and the Laplacian of F(rc)
10.1021/jp002113a CCC: $19.00 © 2000 American Chemical Society Published on Web 10/24/2000
Model Structures for the Silica Polymorph Coesite are found by evaluating the eigenvalues and eigenvectors of the Hessian matrix of F(rc), Hi,j ) ∂2F(rc)/∂xi∂xj, (i,j ) 1, 3). The three eigenvalues of Hi,j are denoted λi (i ) 1, 3), where the Laplacian of F(rc), ∇ 2F(rc) ) λ1 + λ2 + λ3. Both |λ1| and |λ2| define the curvatures of F(r) at rc measured perpendicular to the bond path, and λ3 defines the curvature of F(rc) measured parallel to the bond path. As such, λ1 and λ2 measure the local concentration of F(rc) perpendicular to the bond path and λ3 measures the local depletion of F(rc) parallel to the path in the directions of the nuclei of the bonded atoms. As ∇2F(rc) is equal to the sum of the three eigenvalues, it provides a measure1 of the total local concentration or the local depletion of electron density at rc. When ∇2F(rc) > 0, then F(r) is said to be locally depleted at rc, but when ∇2F(rc) < 0, then F(r) is said to be locally concentrated at rc. However, if for some reason ∇2F(rc) is positive and if it increases in value with decreasing bond length such that the value of the electron density is locally depleted at rc, it does not necessarily follow that F(rc) will decrease in value. An electron density distribution formed by superimposing a set of spherically averaged noninteracting, ground state, static, neutral atoms each clamped at the position that they occupy in a crystal (molecule) is referred to as a procrystal (promolecule) representation of the electron density distribution. In general, the F(rc) values and the curvatures of F(rc) of such a distribution are typically smaller than those observed. However, the procrystal (promolecule) radii of the atoms defined by a procrystal (promolecule) representation are strikingly similar to the bonded radii defined by an observed or a calculated electron density distribution for the crystal (molecule). In this study, the framework structure and the bcp properties of the electron density distribution of the high pressure silica polymorph coesite were modeled at pressure. In earlier studies,2-4 the bcp properties of the electron density distributions1 were calculated for representative Si2O7 and Si5O16 chargesbalanced moieties of the structure. Not only did the resulting properties agree within ∼10%, on average, with those observed at ambient conditions but also with those calculated for the bulk crystal, using the observed cell dimensions and atomic coordinates. The five nonequivalent SiOSi angles observed for the crystal correlate with the observed SiO bond lengths with the shorter bonds involving the wider angles.5 As the bond lengths decrease, the values of the electron densities, F(rc), at the bcps, rc, calculated for the moieties and the bulk crystal,4 increase slightly in agreement with a population analysis completed in earlier MO calculations for the moieties.5 Also, a virial partitioning of the electron density distribution calculated for the bulk crystal shows that the magnitudes of the net charges on the Si and O ions decrease slightly with increasing angle and decreasing bond length as found in a population analysis.5 First-principles calculations have also been completed for the silica polymorph low quartz at zero pressure and pressures up to ∼26 GPa.6 Not only do the bond lengths and angles generated in the calculations match those observed at ambient conditions to within ∼1%, they also both decrease with the compression of the structure.7 The F(rc) values for the two nonequivalent SiO bonds of the model structures increase slightly as the bonded radius, rb(O), and the net charge on the oxide anion both decrease. It was also observed that the length of the shorter bond decreases at a faster rate than the longer one. With compression, the oxide anions of the structure adopt an arrangement that roughly approximates a close-packed array at high pressures.8 In particular, at 26.4 GPa, the model structure6 adopts an SiOSi angle of 120.4°, only ∼10° wider than would
J. Phys. Chem. B, Vol. 104, No. 45, 2000 10535 TABLE 1: Cell Dimensions at Pressure P (GPa)
a (Å)
b (Å)
c (Å)
β (°)
0.0 0.139 1.695 3.302 6.147 8.520 12.622 17.434
7.099 7.094 7.037 6.979 6.894 6.825 6.733 6.648
12.349 12.340 12.299 12.249 12.180 12.118 12.021 11.919
7.170 7.168 7.145 7.120 7.080 7.049 6.998 6.933
120.420 120.422 120.639 120.694 120.938 120.968 121.145 121.262
be adopted by an ideal close-packed structure of oxide anions in which one-quarter of the available tetrahedral voids are occupied by Si cations.8 At ambient conditions and at higher pressures, bond path ridges of electron density and bond critical points were found between a few of the intertetrahedral oxide anions of the model structure, features that have been ascribed to bonded interactions between the oxide anions.9,10 Despite their closer contacts, none was found between the intratetrahedral oxide anions. To assess how the structure and the bcp properties of the electron density distribution of coesite are predicted to behave under compression, first-principles calculation were completed in this study for pressures in excess of 17 GPa. Inasmuch as the coesite structure has been determined at pressures up to ∼5 GPa,11 it will be of interest to see, as calculated for quartz, whether the SiO bond lengths and the SiOSi angles of the model structure both decrease with compression, whether the shorter bonds and the narrower angles are more compressible than the longer bonds and the wider angles, whether the silicate tetrahedra remain largely undistorted, and whether the SiO bond lengths and SiOSi angles decrease with increasing pressure, with the bulk of the compression being accommodated by relatively large changes in the variable SiOSi angles. It will be of interest to see how the topological properties of the electron density distribution respond to pressure and whether bond path ridges of electron density are displayed between the oxide anions as reported for the low quartz model structure6 and other oxides9,10,12 and whether they persist or increase in number with compression. Finally, the ridges will be examined to see whether they can be related to local features of the geometry of the structure. Model Structure Calculations Because of the experimental problems encountered in measuring accurate X-ray diffraction data at high pressures, the observed electron density distributions of coesite and their bcp properties1 have yet to be measured at high pressures. Indeed, given the difficulty of obtaining accurate data at ambient conditions, it is questionable whether accurate electron density distributions and reliable second derivative properties of the distributions will be obtained at high pressures in the near future. In an effort to fulfill this need, first-principles total-energy calculations were undertaken to obtain a set of geometry optimized model structures for coesite determined for a prescribed set of unit cell volumes. In addition, the electron density distributions were calculated, and their bcp properties were determined. In the calculations, each structure was fully geometrically optimized, assuming C2/c space group symmetry and using the Vienna ab initio simulation program (VASP) developed at the Technische Universita¨t Wien, Austria, by Kresse and Furthmu¨ller13 (more details given elsewhere6). The resulting unit cell dimensions and the fractional coordinates of the Si and O atoms in the asymmetric unit of the model structures, calculated for pressures up to 17.4 GPa, are given in Tables 1
10536 J. Phys. Chem. B, Vol. 104, No. 45, 2000
Gibbs et al.
TABLE 2: Fractional Coordinates at Pressure
TABLE 3: SiO Bond Lengths at Pressure
atom
P (GPa)
x
y
z
Si1 Si2 O2 O3 O4 O5 Si1 Si2 O2 O3 O4 O5 Si1 Si2 O2 O3 O4 O5 Si1 Si2 O2 O3 O4 O5 Si1 Si2 O2 O3 O4 O5 Si1 Si2 O2 O3 O4 O5 Si1 Si2 O2 O3 O4 O5 Si1 Si2 O2 O3 O4 O5
0.0
0.13885 0.50709 0.50000 0.73754 0.31326 0.47946 0.13876 0.50701 0.50000 0.73761 0.31313 0.47921 0.13718 0.50834 0.50000 0.74429 0.31764 0.47533 0.13586 0.50821 0.50000 0.74758 0.31866 0.47151 0.13388 0.50918 0.50000 0.75494 0.32306 0.46628 0.13239 0.50861 0.50000 0.75808 0.32402 0.46212 0.13082 0.50879 0.50000 0.76432 0.32767 0.45672 0.12905 0.50873 0.50000 0.76982 0.33116 0.45164
0.10900 0.15749 0.11408 0.12545 0.10282 0.28798 0.10904 0.15739 0.11390 0.12545 0.10274 0.28795 0.10994 0.15714 0.11092 0.12889 0.10169 0.28778 0.11054 0.15687 0.10870 0.13079 0.10081 0.28775 0.11166 0.15657 0.10530 0.13427 0.09979 0.28758 0.11245 0.15619 0.10287 0.13644 0.09896 0.28741 0.11346 0.15583 0.09931 0.13943 0.09796 0.28741 0.11460 0.15550 0.09584 0.14244 0.09722 0.28738
0.07178 0.54198 0.75000 0.56302 0.32719 0.52475 0.07172 0.54201 0.75000 0.56313 0.32715 0.52458 0.07041 0.54394 0.75000 0.56957 0.32572 0.52837 0.06987 0.54465 0.75000 0.57247 0.32524 0.52958 0.06800 0.54659 0.75000 0.57987 0.32347 0.53234 0.06715 0.54728 0.75000 0.58285 0.32277 0.53341 0.06536 0.54885 0.75000 0.58971 0.32127 0.53536 0.06298 0.54999 0.75000 0.59635 0.31925 0.53660
0.139
1.695
3.302
6.147
8.520
12.622
17.434
and 2, respectively. The a-cell edge of the model structure at zero pressure is 0.5% smaller than that observed at ambient conditions, the b-cell edge is 0.15% smaller, the c-cell edge is 0.05% smaller, and the β angle is 0.05% larger, resulting in a unit cell volume that is 1.2% smaller than observed. The SiO bond lengths of the model structures are given in Table 3 and the SiOSi angles are given in Table 4. The total energy, E, was computed as a function of unit cell volume, V, using the generalized gradient approximation to the quantum mechanical energy functional. The equation of state parameters were obtained by fitting a Birch-Murnaghan expression for the energy to the E-V points. The resulting set of total energies together with the set of unit cell volumes were used to compute the zero pressure isothermal bulk modulus, 87.3 GPa, and its pressure derivative, 4.5, as compared with the observed values of 99.8 (( 12.4 GPa) and 6.3 (( 1.2), respectively.14 In the calculations, the Si1O1Si1 angle remained straight over the pressure range studied, and there was no evidence to indicate a change is space group symmetry despite the assumption in several high-pressure calculations of P1 space group symmetry. The CRYSTAL95 program15 was used to generate analytical
P (GPa) R(Si1O1) R(Si1O3) R(Si1O4) R(Si1O5) R(Si2O2) R(Si2O3) R(Si2O4) R(Si2O5) R(Si1O1) R(Si1O3) R(Si1O4) R(Si1O5) R(Si2O2) R(Si2O3) R(Si2O4) R(Si2O5) R(Si1O1) R(Si1O3) R(Si1O4) R(Si1O5) R(Si2O2) R(Si2O3) R(Si2O4) R(Si2O5) R(Si1O1) R(Si1O3) R(Si1O4) R(Si1O5) R(Si2O2) R(Si2O3) R(Si2O4) R(Si2O5)
0.0
0.139
1.695
3.302
R(SiO) (Å)
P (GPa)
R(SiO) (Å)
1.5920 1.6138 1.6113 1.6205 1.6095 1.6142 1.6047 1.6205 1.5910 1.6134 1.6108 1.6198 1.6088 1.6136 1.6041 1.6203 1.5869 1.6126 1.6092 1.6190 1.6070 1.6128 1.6014 1.6191 1.5806 1.6100 1.6061 1.6168 1.6046 1.6107 1.5976 1.6184
6.147
1.5737 1.6076 1.6024 1.6140 1.6005 1.6092 1.5925 1.6165 1.5675 1.6046 1.5996 1.6110 1.5967 1.6062 1.5877 1.6146 1.5584 1.6019 1.5948 1.6074 1.5915 1.6037 1.5807 1.6123 1.5504 1.5979 1.5897 1.6017 1.5860 1.6004 1.5733 1.6084
8.520
12.622
17.434
TABLE 4: SiOSi Angles at Pressure P (GPa)
