- Email: [email protected]
Effect of basis set and of electronic correlation on ab initio calculations on silica rings
Journal of Non-Crystalline Solids 277 (2000) 106±113
www.elsevier.com/locate/jnoncrysol
Eect of basis set and of electronic correlation on ab initio calculations on silica rings J.M. Nedelec a,b,*, L.L. Hench b a
Laboratoire des Mat eriaux Inorganiques, Universit e Blaise Pascal, ESA 6002, 24 Avenue des Landais, F 63 177 Aubi ere cedex, France b Department of Materials, Imperial College of Science Technology and Medicine, Prince Consort Road, London SW7 2BP, UK Received 17 May 2000
Abstract The eect of the choice of basis set and of electronic correlation on ab initio calculations on silica rings have been investigated. Silica rings with size varying between 2 and 6 have been modeled with respect to geometry and energetics. It appears from these results that the 6-31G* basis set is required in order to obtain precis results and that any further improvement of this set is not necessary. Electronic correlation contribution has been estimated thank to DFT calculations. It is concluded that this contribution is very important and cannot be neglected. These two main results build the basis for further modeling of silica based systems including the study of the densi®cation of silica gels. Ó 2000 Elsevier Science B.V. All rights reserved. PACS: 3115.Ar; 3115 Ct; 6143.Bn; 7123.Cq
1. Introduction Silica gels have been recognized as an interesting host matrix for many applications. The mechanical, thermal and optical properties of silica make this material best for numerous optical applications. Doping of porous gel-silica with a secondary active phase gives remarkable properties to the doped hybrid materials. The dopant can include organic dyes to produce solid-state lasers [1] and semiconductor nanoparticles for non-linear optics [2,3]. Use of other dopants gives rise to applications in ®elds such as sensors, optical telecommunications, optical storage, solar cells, etc. * Corresponding author. Tel.: +33-4 73 40 76 47; fax: +33-4 73 40 71 08. E-mail address: [email protected] (J.M. Nedelec).
During processing of silica gels matrices, the densi®cation step is crucial. In eect, within this step, the textural properties of the gels are de®ned and stabilized so that the porous matrices are insensitive to moisture and humidity. Furthermore, the precise control of this densi®cation procedure allows the characteristics of the material such as the density, the pore volume and the pore size distribution to be varied. By modifying the densi®cation temperature, time and atmosphere, one can customize the matrices for a speci®c application. Although the importance of densi®cation is well known, there still exists a lack of understanding of densi®cation mechanisms at an atomic or molecular level. Quantum chemistry can, at least partially, solve this problem. In eect models of amorphous silica have been developed [4±8] and densi®cation pathways have been proposed [9±11]. In most cases, calculations were limited either to
0022-3093/00/$ - see front matter Ó 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 3 0 9 3 ( 0 0 ) 0 0 3 0 6 - 9
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
small clusters or to semi-empirical calculations. These calculations usually required a lot of time on high-capacity computers. Recent developments of microprocessors and calculation codes now allow high level ab initio calculations of clusters of molecules to be run on desktop computers. One aim of this work is to demonstrate that quantum chemistry calculations on small molecules with few electrons can be performed at low cost on a desktop personal computer. Recently [12], the authors have shown that ab initio modeling of silica using a cluster approach was ecient. This previous paper was limited to only one basis set and the electronic correlation was neglected. In the present work, we have used the Spartan ProÓ package from Wavefunction Inc., which provides an interface to various semi-empirical and ab initio Hartree±Fock and DFT methods through a convivial WindowsTM environment. We have performed calculations on silica rings of various sizes with dierent methods studying the eect of the choice of the basis set and of electronic correlation (via DFT calculations) on the structure and energetics of these rings. Limitations of our system led us to restrict DFT calculations to 2±4 membered rings. Nevertheless, these small ring structures play a major role in the densi®cation of silica gels as demonstrated previously [13], and their modeling is essential. 2. Theory The common base to quantum chemistry calculations is solving Schr odinger's equation. In order to do so, several simpli®cations have to be made. Detailed information on this aspect can be found in numerous excellent reviews. Let us focus here on the dierent basis sets [14] available for ab initio Hartree±Fock calculations and on the use of density functional theory (DFT) to account for electronic correlation. 2.1. Basis sets For practical reasons, Hartree±Fock calculations use gaussian-type wavefunctions, the most popular being minimal (STO-3G) and split-valence
107
(3-21G, 6-31G and 6-311G) basis sets supplemented or not by sets of polarization and±or diffuse functions. The simplest basis set is the STO-3G set [15], in which each of the basis functions is expanded in terms of three gaussians. This set does not describe correctly atoms with aspherical environments and electron distributions between nuclei because of its atom-centered character. To improve the minimal basis set, it has been proposed to provide two sets of valence basis functions. This kind of split-valence sets represents the inner-shell atomic orbitals (a.o.) by a single set of functions and the valence atomic orbitals by two sets of functions. The splitvalence sets used in this study are the 3-21G and 6-31G sets. Each inner-shell a.o. in the 3-21G set is represented by a single function which is written in terms of three gaussians, while the valence a.o. are written in terms of expansion of two and one gaussians, respectively. 6-31G sets are constructed in the same way with inner-shell orbitals represented in terms of six gaussians and valence orbitals split into three and one gaussian components. One further improvement of minimal or splitvalence basis sets is the introduction of a polarization basis set [16]. The latter provides d-type functions on main group atoms. The simpler polarization basis sets are 3-21G* and 6-31G* and are derived from 3-21G and 6-31G sets, respectively, by adding a set of six d-type gaussians for each heavy atom. A further improvement, called 6-31G**, is identical to 6-31G* except that it provides p-type polarization functions for hydrogen atoms. A ®nal re®nement of the basis set implies the introduction of diuse functions [17] for heavy atoms which are supplemented with diuse s- and p-type functions. It leads to 6-31 G and 6-31 G sets. A higher degree of splitting is possible with the 6-311G basis set producing the resulting 6-311G*, 6-311G** and 6-311 G sets, but they have not been used in this work because of the high cost of these calculations. 2.2. DFT calculations Hartree±Fock models treat the motions of individual electrons as independent of one another.
108
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
This leads to overestimation of the electron±electron repulsion energy. The correlation energy is de®ned as the dierence between the real energy and the Hartree±Fock energy. Several methods have been developed to account for this correlation energy among which the most popular is probably the Mùller±Plesset perturbation technique. Another solution to treat the correlation in many electrons systems is the density functional theory (DFT); see for example [18] for a review. This theory is based on the fact that all the properties of a system (in particular the energy) can be derived from the knowledge of its electronic density. The simplest DFT models are called local spin density and referred to as SVWN (Slater, Vosko, Wilk and Nusair). An improvement can be made by introducing an explicit dependence on the gradient of the electron density, in addition to the density itself. Within Spartan, it leads to the perturbative Becke-Perdew model (pBP) [19,20]. The DFT models available in Spartan use tabulated atomic solutions supplemented by d-type functions on heavy atoms instead of gaussian basis sets. Numerical polarization (*) and full polarization (**) basis sets are referred to as DN* and DN**, respectively. In practice, the methods available are SVWN/DN*, SVWN/DN**, pBP/DN* and pBP/DN**. 3. Results In this work, amorphous silica has been modeled using a cluster approach with ring size varying between 2 and 6 tetrahedra, all tetrahedra being bonded together with O±Si±O siloxane bonds. The rings have been modeled using a Hartree±Fock ab initio method with several basis sets previously described, namely STO-3G, 3-21G*, 6-31G*, 6-31G** and 6-31 G . All rings were terminated with silanol groups leading to the formulation H2n Sin O3n with n varying between 2 and 6. The geometry of the rings has been optimized and the total electronic energy determined. The geometrical parameters (Si±O±Si and O±Si±O bond angles and Si±OH and Si±O±Si bond length) can be found in Table 1 together with the energy per tetrahedral unit for all the basis sets used. Results obtained
with the 3-21G* basis set have been taken from [12]. In order to compare results from the dierent basis sets, we have calculated the average of the geometrical parameters over the dierent ring sizes taking into account the size of the ring and its relative concentration in amorphous silica. The relative concentration of the dierent rings has been taken from the work of Bell and Dean [21]. 3.1. Hartree±Fock calculations ± Eect of the basis set 3.1.1. Bond angles The evolution of the Si±O±Si bond angle for the dierent basis sets as a function of n is presented in Fig. 1. The bond angle shows little dependence on the choice of the basis set for small rings. Larger bond angles dependence is obtained for n 4±6. Results obtained for 6-31G*, 6-31G** and 631 G are very similar and the curves are superimposed. Results for the O±Si±O bond angle are illustrated in Fig. 2 and show very little discrepancy between all the basis sets. Once again the curves obtained for 6-31G*, 6-31G** and 6-31 G are equivalent. 3.1.2. Bond length The evolution of the non-bridging Si±OH bond length is shown in Fig. 3. While results for the three more sophisticated basis sets are still very similar, a very dierent behavior is observed for STO-3G and 3-21G*. The minimal basis set shows a maximum length of Si±OH bonds for the 5-membered ring which is unjusti®ed. The evolution of the bridging Si±O±Si bond length depicted in Fig. 4 is similar to the one obtained for non-bridging oxygen bonds with a maximum for the 5-membered ring for the STO3G set. 3.1.3. Energetics The energy per tetrahedra for the dierent size rings is shown in Fig. 5. Scales have been adjusted to visualize all the curves on the same graph. The values are consequently only representative of the dierences. Precise values are given in Table 1. The trend for all the basis sets is the same. The highly
a
2 3 4 5 6 B&D
2 3 4 5 6 B&D
2 3 4 5 6 B&D
2 3 4 5 6
O±Si±O
Si±(OH)
Si±(OSi)
Energy
)508.383 )508.413 )508.418 )508.41 )508.422
1.676 1.619 1.618 1.644 1.619 1.625
1.648 1.649 1.649 1.675 1.648 1.655
87.17 107.02 107.43 112.46 108.80 109.54
92.83 131.03 132.23 140.71 132.08 134.33
STO-3G
)512.25 )512.281 )512.294 )512.296 )512.303
1.61 1.61 1.62 1.62 1.62 1.620
1.66 1.63 1.62 1.61 1.62 1.618
88 103.5 107 109 107 107.47
92 136 142 157 155 153.46
3-21G*
)514.857 )514.882 )514.887 )514.888 )514.89
1.659 1.628 1.625 1.621 1.619 1.621
1.616 1.623 1.624 1.627 1.627 1.627
88.21 105.85 108.49 111.88 109.23 109.77
91.79 133.95 138.91 140.48 145.23 142.96
6-31G*
)514.872 )514.897 )514.901 )514.903 )514.904
1.659 1.628 1.625 1.625 1.619 1.621
1.614 1.620 1.621 1.621 1.625 1.623
88.21 105.88 108.57 111.57 109.24 109.70
91.79 133.93 139.17 140.80 145.61 143.31
6-31G**
)514.867 )514.893 )514.897 )514.898 )514.902
1.659 1.627 1.625 1.625 1.619 1.621
1.612 1.623 1.621 1.621 1.625 1.623
88.10 105.84 108.58 111.57 109.24 109.70
91.91 133.97 139.17 140.80 145.61 143.31
6-31 + G
1.660 1.628 1.630
1.619 1.620 1.624
90.27 106.88 109.32
)513.948 )513.975 )513.979 ) )
) ) )
) ) )
) ) )
) ) )
89.73 133.06 128.62
SVWN/ DN*
1.660 1.629 1.631
)513.957 )513.986 )513.99 ) )
) ) )
) ) )
1.618 1.619 1.622
90.25 106.88 109.26 ) ) )
89.75 133.07 128.85 ) ) )
SVWN/ DN**
1.700 1.650 1.648
1.685 1.641 1.642
84.32 107.16 110.07
)516.708 )516.732 )516.732 ) )
) ) )
) ) )
) ) )
) ) )
95.68 132.79 133.97
PBP/DN*
1.682 1.650 1.649
1.639 1.639 1.641
90.29 107.15 109.99
)516.718 )516.743 )516.744 ) )
) ) )
) ) )
) ) )
) ) )
89.71 132.79 134.19
PBP/ DN**
B&D stands for the average value obtained by weighing the values obtained for the rings by their concentration derived from the Bell and Dean distribution.
2 3 4 5 6 B&D
Si±O±Si
n
Table 1 Geometrical parameters and energies per tetrahedral unit for the various rings calculated with the dierent basis setsa
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113 109
110
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
Fig. 1. Evolution of the Si±O±Si bond angle as a function of the ring size n, for the dierent basis sets.
Fig. 2. Evolution of the O±Si±O bond angle as a function of the ring size n, for the dierent basis sets.
Fig. 3. Evolution of the Si±OH bond length as a function of the ring size n, for the dierent basis sets.
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
111
Fig. 4. Evolution of the Si±OSi bond length as a function of the ring size n, for the dierent basis sets.
Fig. 5. Energy per tetrahedral unit for the dierent rings modeled with the dierent basis sets.
stressed 2-membered ring has a very high energy and the energy decreases with ring size, the bigger ring being the more stable. 3.2. DFT calculations ± correlation energy DFT calculations have been performed on rings with 2±4 tetrahedral units with four dierent basis sets: SVWN/DN*, SVWN/DN**, pBP/DN* and pBP/DN**. All the results are shown in Table 1. The geometries are quite dierent from the Hartree±Fock models especially for the perturbative Becke Perdew models. As observed for the non-correlated calculations, the 2-membered ring is unstable with the highest energy, while the biggest ring ( n 4
for the DFT calculations) is the more stable. More interesting is the comparison of the values themselves. The energies from the local DFT models appear to be between the energies of the 3-21G and the 6-31G models. But for the more sophisticated DFT models (pBP) the energies are much lower by some 2 a.u. (1255 kCal/mol !!) lower than the 6-31G models.
4. Discussion 4.1. Eect of the basis set The geometrical parameters derived from the calculations have to be compared with the
112
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
experimental values obtained for amorphous silica. We have used the results from Mozzi and Warren [22] who have analyzed X-ray diraction pair function of condensed phase amorphous silica. These authors found a Si±O±Si bond length and a most probable value of Si±O±Si of 1.62 A bond angle of 144°. The value of the O±Si±O angle depends very little on the basis set and from the 6-31G* set it reaches a value corresponding to the ideal value of 109:45°, corresponding to a perfect tetrahedral SiO4 unit. These results indicate, as stated before, that the basic tetrahedral unit is not distorted and that the stress resulting from the ring formation is accommodated by varying the inter-tetrahedral angle Si±O±Si. The value obtained for the Si±O±Si angle is very good for the 6-31G* basis set and the more sophisticated sets. It is around 143° compared to 144° for the experimental value. The value for the Si±O±Si bond length is also in very good agree except ment with the experimental result (1:62 A) for the STO3G set. In general, the quality of the simulation is increased from STO3G to 3-21G to 6-31G. STO3G is clearly not satisfying and 6-31G* gives very good results. Any further improvement of the 6-31G* basis set is not necessary since the changes are insigni®cant. 4.2. Correlation energy The electronic correlation energy contribution is important even if the general trend is the same as obtained without taking it into account (the bigger the ring, the more stable it is). Nevertheless, for precise energetic calculations, this energy cannot be neglected as a discrepancy of as high as 2 a.u. can be observed. The geometry derived from the DFT calculations is not very good with comparison to experimental data. One solution would be to develop customized basis sets for DFT modeling of silica rings. The HF calculations with 6-31G* basis set could serve as a basis for such improvement. As stated in the introduction of this work, the goal of these calculations is to obtain models of silica structural units that can be subsequently rearranged to simulate densi®cation of
silica gels during thermal treatment. To achieve this objective, precise energetics calculations are required. It follows from this work that the electronic correlation energy has to be taken into account. The DFT method appear to be the simplest way to do so even if it requires improvements. 5. Conclusion A thorough analysis of the eects of basis set and of correlation energy on the ab initio calculations on silica rings has been performed. It has been shown that HF calculations require the 6-31G* level to obtain accurate geometries and that any further improvement beyond this basis set is not useful. Furthermore, evaluation of the correlation energy through DFT calculations led to the conclusion that this energy cannot be neglected and must be taken into account for precise energetics calculations on silica systems. These two main conclusions provide the basis for further modeling of silica gels, in particular, the study of the evolution of the microstructure of these materials with respect to thermal treatment. Acknowledgements The authors gratefully acknowledge the support of the United Kingdom Engineering and Physical Sciences Research Council under projects 86265 and 78185. References [1] R. Reisfeld, J. Lumin. 72±74 (1997) 7. [2] J. Butty, Y.Z. Hu, N. Peyghambarian, Y.H. Kao, J.D. Mackenzie, Appl. Phys. Lett. 67 (1995) 2672. [3] B. Capoen, T. Gacoin, J.M. Nedelec, S. Turrell, M. Bouazaoui, J. Mater. Sci. (2000) in press. [4] D.B. Boyd, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, VCH, New York, 1990, p. 321. [5] B.W.H. van Beest, J. Verbeek, R.A. Van Santen, Catal. Lett. 1 (1988) 147. [6] A.C. Hess, P.F. McMillan, M. O'Keefe, J. Chem. Phys. 90 (1986) 5661. [7] J.K. West, L.L. Hench, J. Non-Cryst. Solids 180 (1994) 11. [8] J.A. Harkless, D.K. Stillinger, F.H. Stillinger, J. Phys. Chem 100 (1996) 1098.
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113 [9] [10] [11] [12]
J.K. West, L.L. Hench, J. Mater. Sci. 29 (1994) 3601. J.K. West, L.L. Hench, J. Mater. Sci. 29 (1994) 5808. J.K. West, L.L. Hench, J. Mater. Sci. 30 (1995) 6281. J.M. Nedelec, L.L. Hench, J. Non-Cryst. Solids 255 (1999) 163. [13] J.M. Nedelec, M. Bouazaoui, S. Turrell, J. Non-Cryst. Solids 243 (1999) 209. [14] D. Feller, E.R. Davidson, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, vol. 1, VCH, New York, 1990, p. 1. [15] W.J. Hehre, R.F. Stewart, J.A. Pople, J. Chem. Phys. 51 (1969) 2657.
113
[16] P.C. Hariharan, J.A. Pople, Chem. Phys. Lett. 66 (1972) 217. [17] T. Clark, J. Chandrasekhar, K. Spitznagel, P.v.R. Schleyer, J. Comput. Chem. 4 (1983) 294. [18] R.G. Parr, W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University, Oxford, 1989. [19] A.D. Becke, Phys. Rev. A 38 (1988) 3089. [20] J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [21] R.J. Bell, P. Dean, Philos. Mag. 25 (1972) 1381. [22] R.L. Mozzi, B.E. Warren, J. Appl. Crystallogr. 2 (1969) 164.
www.elsevier.com/locate/jnoncrysol
Eect of basis set and of electronic correlation on ab initio calculations on silica rings J.M. Nedelec a,b,*, L.L. Hench b a
Laboratoire des Mat eriaux Inorganiques, Universit e Blaise Pascal, ESA 6002, 24 Avenue des Landais, F 63 177 Aubi ere cedex, France b Department of Materials, Imperial College of Science Technology and Medicine, Prince Consort Road, London SW7 2BP, UK Received 17 May 2000
Abstract The eect of the choice of basis set and of electronic correlation on ab initio calculations on silica rings have been investigated. Silica rings with size varying between 2 and 6 have been modeled with respect to geometry and energetics. It appears from these results that the 6-31G* basis set is required in order to obtain precis results and that any further improvement of this set is not necessary. Electronic correlation contribution has been estimated thank to DFT calculations. It is concluded that this contribution is very important and cannot be neglected. These two main results build the basis for further modeling of silica based systems including the study of the densi®cation of silica gels. Ó 2000 Elsevier Science B.V. All rights reserved. PACS: 3115.Ar; 3115 Ct; 6143.Bn; 7123.Cq
1. Introduction Silica gels have been recognized as an interesting host matrix for many applications. The mechanical, thermal and optical properties of silica make this material best for numerous optical applications. Doping of porous gel-silica with a secondary active phase gives remarkable properties to the doped hybrid materials. The dopant can include organic dyes to produce solid-state lasers [1] and semiconductor nanoparticles for non-linear optics [2,3]. Use of other dopants gives rise to applications in ®elds such as sensors, optical telecommunications, optical storage, solar cells, etc. * Corresponding author. Tel.: +33-4 73 40 76 47; fax: +33-4 73 40 71 08. E-mail address: [email protected] (J.M. Nedelec).
During processing of silica gels matrices, the densi®cation step is crucial. In eect, within this step, the textural properties of the gels are de®ned and stabilized so that the porous matrices are insensitive to moisture and humidity. Furthermore, the precise control of this densi®cation procedure allows the characteristics of the material such as the density, the pore volume and the pore size distribution to be varied. By modifying the densi®cation temperature, time and atmosphere, one can customize the matrices for a speci®c application. Although the importance of densi®cation is well known, there still exists a lack of understanding of densi®cation mechanisms at an atomic or molecular level. Quantum chemistry can, at least partially, solve this problem. In eect models of amorphous silica have been developed [4±8] and densi®cation pathways have been proposed [9±11]. In most cases, calculations were limited either to
0022-3093/00/$ - see front matter Ó 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 3 0 9 3 ( 0 0 ) 0 0 3 0 6 - 9
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
small clusters or to semi-empirical calculations. These calculations usually required a lot of time on high-capacity computers. Recent developments of microprocessors and calculation codes now allow high level ab initio calculations of clusters of molecules to be run on desktop computers. One aim of this work is to demonstrate that quantum chemistry calculations on small molecules with few electrons can be performed at low cost on a desktop personal computer. Recently [12], the authors have shown that ab initio modeling of silica using a cluster approach was ecient. This previous paper was limited to only one basis set and the electronic correlation was neglected. In the present work, we have used the Spartan ProÓ package from Wavefunction Inc., which provides an interface to various semi-empirical and ab initio Hartree±Fock and DFT methods through a convivial WindowsTM environment. We have performed calculations on silica rings of various sizes with dierent methods studying the eect of the choice of the basis set and of electronic correlation (via DFT calculations) on the structure and energetics of these rings. Limitations of our system led us to restrict DFT calculations to 2±4 membered rings. Nevertheless, these small ring structures play a major role in the densi®cation of silica gels as demonstrated previously [13], and their modeling is essential. 2. Theory The common base to quantum chemistry calculations is solving Schr odinger's equation. In order to do so, several simpli®cations have to be made. Detailed information on this aspect can be found in numerous excellent reviews. Let us focus here on the dierent basis sets [14] available for ab initio Hartree±Fock calculations and on the use of density functional theory (DFT) to account for electronic correlation. 2.1. Basis sets For practical reasons, Hartree±Fock calculations use gaussian-type wavefunctions, the most popular being minimal (STO-3G) and split-valence
107
(3-21G, 6-31G and 6-311G) basis sets supplemented or not by sets of polarization and±or diffuse functions. The simplest basis set is the STO-3G set [15], in which each of the basis functions is expanded in terms of three gaussians. This set does not describe correctly atoms with aspherical environments and electron distributions between nuclei because of its atom-centered character. To improve the minimal basis set, it has been proposed to provide two sets of valence basis functions. This kind of split-valence sets represents the inner-shell atomic orbitals (a.o.) by a single set of functions and the valence atomic orbitals by two sets of functions. The splitvalence sets used in this study are the 3-21G and 6-31G sets. Each inner-shell a.o. in the 3-21G set is represented by a single function which is written in terms of three gaussians, while the valence a.o. are written in terms of expansion of two and one gaussians, respectively. 6-31G sets are constructed in the same way with inner-shell orbitals represented in terms of six gaussians and valence orbitals split into three and one gaussian components. One further improvement of minimal or splitvalence basis sets is the introduction of a polarization basis set [16]. The latter provides d-type functions on main group atoms. The simpler polarization basis sets are 3-21G* and 6-31G* and are derived from 3-21G and 6-31G sets, respectively, by adding a set of six d-type gaussians for each heavy atom. A further improvement, called 6-31G**, is identical to 6-31G* except that it provides p-type polarization functions for hydrogen atoms. A ®nal re®nement of the basis set implies the introduction of diuse functions [17] for heavy atoms which are supplemented with diuse s- and p-type functions. It leads to 6-31 G and 6-31 G sets. A higher degree of splitting is possible with the 6-311G basis set producing the resulting 6-311G*, 6-311G** and 6-311 G sets, but they have not been used in this work because of the high cost of these calculations. 2.2. DFT calculations Hartree±Fock models treat the motions of individual electrons as independent of one another.
108
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
This leads to overestimation of the electron±electron repulsion energy. The correlation energy is de®ned as the dierence between the real energy and the Hartree±Fock energy. Several methods have been developed to account for this correlation energy among which the most popular is probably the Mùller±Plesset perturbation technique. Another solution to treat the correlation in many electrons systems is the density functional theory (DFT); see for example [18] for a review. This theory is based on the fact that all the properties of a system (in particular the energy) can be derived from the knowledge of its electronic density. The simplest DFT models are called local spin density and referred to as SVWN (Slater, Vosko, Wilk and Nusair). An improvement can be made by introducing an explicit dependence on the gradient of the electron density, in addition to the density itself. Within Spartan, it leads to the perturbative Becke-Perdew model (pBP) [19,20]. The DFT models available in Spartan use tabulated atomic solutions supplemented by d-type functions on heavy atoms instead of gaussian basis sets. Numerical polarization (*) and full polarization (**) basis sets are referred to as DN* and DN**, respectively. In practice, the methods available are SVWN/DN*, SVWN/DN**, pBP/DN* and pBP/DN**. 3. Results In this work, amorphous silica has been modeled using a cluster approach with ring size varying between 2 and 6 tetrahedra, all tetrahedra being bonded together with O±Si±O siloxane bonds. The rings have been modeled using a Hartree±Fock ab initio method with several basis sets previously described, namely STO-3G, 3-21G*, 6-31G*, 6-31G** and 6-31 G . All rings were terminated with silanol groups leading to the formulation H2n Sin O3n with n varying between 2 and 6. The geometry of the rings has been optimized and the total electronic energy determined. The geometrical parameters (Si±O±Si and O±Si±O bond angles and Si±OH and Si±O±Si bond length) can be found in Table 1 together with the energy per tetrahedral unit for all the basis sets used. Results obtained
with the 3-21G* basis set have been taken from [12]. In order to compare results from the dierent basis sets, we have calculated the average of the geometrical parameters over the dierent ring sizes taking into account the size of the ring and its relative concentration in amorphous silica. The relative concentration of the dierent rings has been taken from the work of Bell and Dean [21]. 3.1. Hartree±Fock calculations ± Eect of the basis set 3.1.1. Bond angles The evolution of the Si±O±Si bond angle for the dierent basis sets as a function of n is presented in Fig. 1. The bond angle shows little dependence on the choice of the basis set for small rings. Larger bond angles dependence is obtained for n 4±6. Results obtained for 6-31G*, 6-31G** and 631 G are very similar and the curves are superimposed. Results for the O±Si±O bond angle are illustrated in Fig. 2 and show very little discrepancy between all the basis sets. Once again the curves obtained for 6-31G*, 6-31G** and 6-31 G are equivalent. 3.1.2. Bond length The evolution of the non-bridging Si±OH bond length is shown in Fig. 3. While results for the three more sophisticated basis sets are still very similar, a very dierent behavior is observed for STO-3G and 3-21G*. The minimal basis set shows a maximum length of Si±OH bonds for the 5-membered ring which is unjusti®ed. The evolution of the bridging Si±O±Si bond length depicted in Fig. 4 is similar to the one obtained for non-bridging oxygen bonds with a maximum for the 5-membered ring for the STO3G set. 3.1.3. Energetics The energy per tetrahedra for the dierent size rings is shown in Fig. 5. Scales have been adjusted to visualize all the curves on the same graph. The values are consequently only representative of the dierences. Precise values are given in Table 1. The trend for all the basis sets is the same. The highly
a
2 3 4 5 6 B&D
2 3 4 5 6 B&D
2 3 4 5 6 B&D
2 3 4 5 6
O±Si±O
Si±(OH)
Si±(OSi)
Energy
)508.383 )508.413 )508.418 )508.41 )508.422
1.676 1.619 1.618 1.644 1.619 1.625
1.648 1.649 1.649 1.675 1.648 1.655
87.17 107.02 107.43 112.46 108.80 109.54
92.83 131.03 132.23 140.71 132.08 134.33
STO-3G
)512.25 )512.281 )512.294 )512.296 )512.303
1.61 1.61 1.62 1.62 1.62 1.620
1.66 1.63 1.62 1.61 1.62 1.618
88 103.5 107 109 107 107.47
92 136 142 157 155 153.46
3-21G*
)514.857 )514.882 )514.887 )514.888 )514.89
1.659 1.628 1.625 1.621 1.619 1.621
1.616 1.623 1.624 1.627 1.627 1.627
88.21 105.85 108.49 111.88 109.23 109.77
91.79 133.95 138.91 140.48 145.23 142.96
6-31G*
)514.872 )514.897 )514.901 )514.903 )514.904
1.659 1.628 1.625 1.625 1.619 1.621
1.614 1.620 1.621 1.621 1.625 1.623
88.21 105.88 108.57 111.57 109.24 109.70
91.79 133.93 139.17 140.80 145.61 143.31
6-31G**
)514.867 )514.893 )514.897 )514.898 )514.902
1.659 1.627 1.625 1.625 1.619 1.621
1.612 1.623 1.621 1.621 1.625 1.623
88.10 105.84 108.58 111.57 109.24 109.70
91.91 133.97 139.17 140.80 145.61 143.31
6-31 + G
1.660 1.628 1.630
1.619 1.620 1.624
90.27 106.88 109.32
)513.948 )513.975 )513.979 ) )
) ) )
) ) )
) ) )
) ) )
89.73 133.06 128.62
SVWN/ DN*
1.660 1.629 1.631
)513.957 )513.986 )513.99 ) )
) ) )
) ) )
1.618 1.619 1.622
90.25 106.88 109.26 ) ) )
89.75 133.07 128.85 ) ) )
SVWN/ DN**
1.700 1.650 1.648
1.685 1.641 1.642
84.32 107.16 110.07
)516.708 )516.732 )516.732 ) )
) ) )
) ) )
) ) )
) ) )
95.68 132.79 133.97
PBP/DN*
1.682 1.650 1.649
1.639 1.639 1.641
90.29 107.15 109.99
)516.718 )516.743 )516.744 ) )
) ) )
) ) )
) ) )
) ) )
89.71 132.79 134.19
PBP/ DN**
B&D stands for the average value obtained by weighing the values obtained for the rings by their concentration derived from the Bell and Dean distribution.
2 3 4 5 6 B&D
Si±O±Si
n
Table 1 Geometrical parameters and energies per tetrahedral unit for the various rings calculated with the dierent basis setsa
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113 109
110
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
Fig. 1. Evolution of the Si±O±Si bond angle as a function of the ring size n, for the dierent basis sets.
Fig. 2. Evolution of the O±Si±O bond angle as a function of the ring size n, for the dierent basis sets.
Fig. 3. Evolution of the Si±OH bond length as a function of the ring size n, for the dierent basis sets.
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
111
Fig. 4. Evolution of the Si±OSi bond length as a function of the ring size n, for the dierent basis sets.
Fig. 5. Energy per tetrahedral unit for the dierent rings modeled with the dierent basis sets.
stressed 2-membered ring has a very high energy and the energy decreases with ring size, the bigger ring being the more stable. 3.2. DFT calculations ± correlation energy DFT calculations have been performed on rings with 2±4 tetrahedral units with four dierent basis sets: SVWN/DN*, SVWN/DN**, pBP/DN* and pBP/DN**. All the results are shown in Table 1. The geometries are quite dierent from the Hartree±Fock models especially for the perturbative Becke Perdew models. As observed for the non-correlated calculations, the 2-membered ring is unstable with the highest energy, while the biggest ring ( n 4
for the DFT calculations) is the more stable. More interesting is the comparison of the values themselves. The energies from the local DFT models appear to be between the energies of the 3-21G and the 6-31G models. But for the more sophisticated DFT models (pBP) the energies are much lower by some 2 a.u. (1255 kCal/mol !!) lower than the 6-31G models.
4. Discussion 4.1. Eect of the basis set The geometrical parameters derived from the calculations have to be compared with the
112
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113
experimental values obtained for amorphous silica. We have used the results from Mozzi and Warren [22] who have analyzed X-ray diraction pair function of condensed phase amorphous silica. These authors found a Si±O±Si bond length and a most probable value of Si±O±Si of 1.62 A bond angle of 144°. The value of the O±Si±O angle depends very little on the basis set and from the 6-31G* set it reaches a value corresponding to the ideal value of 109:45°, corresponding to a perfect tetrahedral SiO4 unit. These results indicate, as stated before, that the basic tetrahedral unit is not distorted and that the stress resulting from the ring formation is accommodated by varying the inter-tetrahedral angle Si±O±Si. The value obtained for the Si±O±Si angle is very good for the 6-31G* basis set and the more sophisticated sets. It is around 143° compared to 144° for the experimental value. The value for the Si±O±Si bond length is also in very good agree except ment with the experimental result (1:62 A) for the STO3G set. In general, the quality of the simulation is increased from STO3G to 3-21G to 6-31G. STO3G is clearly not satisfying and 6-31G* gives very good results. Any further improvement of the 6-31G* basis set is not necessary since the changes are insigni®cant. 4.2. Correlation energy The electronic correlation energy contribution is important even if the general trend is the same as obtained without taking it into account (the bigger the ring, the more stable it is). Nevertheless, for precise energetic calculations, this energy cannot be neglected as a discrepancy of as high as 2 a.u. can be observed. The geometry derived from the DFT calculations is not very good with comparison to experimental data. One solution would be to develop customized basis sets for DFT modeling of silica rings. The HF calculations with 6-31G* basis set could serve as a basis for such improvement. As stated in the introduction of this work, the goal of these calculations is to obtain models of silica structural units that can be subsequently rearranged to simulate densi®cation of
silica gels during thermal treatment. To achieve this objective, precise energetics calculations are required. It follows from this work that the electronic correlation energy has to be taken into account. The DFT method appear to be the simplest way to do so even if it requires improvements. 5. Conclusion A thorough analysis of the eects of basis set and of correlation energy on the ab initio calculations on silica rings has been performed. It has been shown that HF calculations require the 6-31G* level to obtain accurate geometries and that any further improvement beyond this basis set is not useful. Furthermore, evaluation of the correlation energy through DFT calculations led to the conclusion that this energy cannot be neglected and must be taken into account for precise energetics calculations on silica systems. These two main conclusions provide the basis for further modeling of silica gels, in particular, the study of the evolution of the microstructure of these materials with respect to thermal treatment. Acknowledgements The authors gratefully acknowledge the support of the United Kingdom Engineering and Physical Sciences Research Council under projects 86265 and 78185. References [1] R. Reisfeld, J. Lumin. 72±74 (1997) 7. [2] J. Butty, Y.Z. Hu, N. Peyghambarian, Y.H. Kao, J.D. Mackenzie, Appl. Phys. Lett. 67 (1995) 2672. [3] B. Capoen, T. Gacoin, J.M. Nedelec, S. Turrell, M. Bouazaoui, J. Mater. Sci. (2000) in press. [4] D.B. Boyd, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, VCH, New York, 1990, p. 321. [5] B.W.H. van Beest, J. Verbeek, R.A. Van Santen, Catal. Lett. 1 (1988) 147. [6] A.C. Hess, P.F. McMillan, M. O'Keefe, J. Chem. Phys. 90 (1986) 5661. [7] J.K. West, L.L. Hench, J. Non-Cryst. Solids 180 (1994) 11. [8] J.A. Harkless, D.K. Stillinger, F.H. Stillinger, J. Phys. Chem 100 (1996) 1098.
J.M. Nedelec, L.L. Hench / Journal of Non-Crystalline Solids 277 (2000) 106±113 [9] [10] [11] [12]
J.K. West, L.L. Hench, J. Mater. Sci. 29 (1994) 3601. J.K. West, L.L. Hench, J. Mater. Sci. 29 (1994) 5808. J.K. West, L.L. Hench, J. Mater. Sci. 30 (1995) 6281. J.M. Nedelec, L.L. Hench, J. Non-Cryst. Solids 255 (1999) 163. [13] J.M. Nedelec, M. Bouazaoui, S. Turrell, J. Non-Cryst. Solids 243 (1999) 209. [14] D. Feller, E.R. Davidson, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, vol. 1, VCH, New York, 1990, p. 1. [15] W.J. Hehre, R.F. Stewart, J.A. Pople, J. Chem. Phys. 51 (1969) 2657.
113
[16] P.C. Hariharan, J.A. Pople, Chem. Phys. Lett. 66 (1972) 217. [17] T. Clark, J. Chandrasekhar, K. Spitznagel, P.v.R. Schleyer, J. Comput. Chem. 4 (1983) 294. [18] R.G. Parr, W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University, Oxford, 1989. [19] A.D. Becke, Phys. Rev. A 38 (1988) 3089. [20] J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [21] R.J. Bell, P. Dean, Philos. Mag. 25 (1972) 1381. [22] R.L. Mozzi, B.E. Warren, J. Appl. Crystallogr. 2 (1969) 164.