Ab initio calculations on borate systems

Proc. Fifth Int. Conf. on Borate Glasses, Crystals and Melts Phys. Chem. Glasses: Eur. J. Glass Sci. Technol. B, August 2006, 47 (4), 441–444

Ab initio calculations on borate systems G. Ferlat,1 L. Cormier, F. Mauri, E. Balan, G. Calas Institut de Mineralogie et de Physique des Milieux Condenses, UMR CNRS-IPGP-Universites Paris 6 et 7 n° 7590, case 115 4 place Jussieu, 75252 Paris cedex 05, France

T. Charpentier Laboratoire Claude Frejacques, URA CNRS-CEA n° 331 CEA Saclay, 91191 Gif-sur-Yvette cedex, France

E. Anglada Departamento de Fisica de la Materia Condensada, Universidad Autonoma de Madrid E-28049 Madrid, Spain This paper is part of an ongoing effort aimed at modelling the structure of liquid and vitreous B2O3 within a first principles framework. Using density-functional theory, pseudoatomic orbitals basis sets and pseudopotentials, we carried out a series of calculations on known borate systems, namely boric acid B(OH)3 and the two polymorphic B2O3 crystals. The results obtained show that rather small basis sizes can compete with the much more expensive planewave ones. This paves the way for an efficient description of the disordered phases using ab-initio molecular dynamics simulations.

1. Introduction Borate glasses, crystals and melts are systems of considerable interest from both fundamental and technological point of views. Boron trioxide B2O3, as a key component of such systems, has attracted much attention. Regarding its structural description, it is however striking that even in such a simple case, there remain much debated issues such as the existence of phase transitions between amorphous,(1,2) and even liquid phases(3) or the importance of medium range superstructural units, i.e. the so-called boroxol (B2O6) rings, below the glass transition.(4,5) Regarding the latter issue, the theoretical description of B2O3 glasses, initiated more than two decades ago by molecular dynamics (MD) simulations,(6–8) appeared to be quite challenging. Indeed, while most of the experimental information tends to establish that a high proportion (~80%) of the boron atoms are involved in boroxol units, the numerical investigations predicted so far a much lower figure (in between 0 and 30%). This failure has been attributed to the use of empirical force fields and possibly to the use of unrealistic quenching rates in the simulated samples. In fact, simulations using potentials restricted to pair interactions did not succeed at all in reproducing boroxos-containing structures and it appeared that the inclusion of higher order (three- and even fourbody) interactions are required in order to get some of the BO3 triangles to connect into boroxol rings.(9–11) It has been suggested that the potential should include forces necessary to account for the ring stabilisation 1

Corresponding author

energy which arises from delocalised π-bonding.(4) There has been several ab-initio total energy calculations undertaken on finite-size clusters(12–14) or crystals(15–18) and such types of calculations have been used to derive empirical potentials for MD simulations.(10,18,19) However, the obtained fraction of boroxols was in any case below 30%. There has not been yet any full ab-initio MD simulations of bulk liquid and glassy B2O3 although it is clear from the above considerations that a first-principles description is highly desirable. This lack stems from the highly time demanding needs of such simulations which until recently have limited the size of the simulated system to several tens of atoms. In addition, the accessible time-scale over which the simulations are carried out are typically of the order of few picoseconds. These limitations in size and time severely hinder the simulation of a realistic glass. However, significant progress toward more realistic system sizes (and simulation times) have been made recently, among which the development of the SIESTA code.(20) This code, based on the the density functional theory (DFT), allows for a very efficient resolution of the Kohn–Sham equations thanks to the use of pseudo-atomic orbitals basis; it is thus very adapted to simulations of disordered systems for which large unit cells are required. In addition, the solution to the eigenvalue problem can be performed either by diagonalisation or with the linear scaling solver of Kim et al.(21) This allows for simulation times scaling essentially like the number of particles instead of the usual cubic behaviour (see for instance Fernández-Serra et al(22) for a successful

Physics and Chemistry of Glasses: European Journal of Glass Science and Technology Part B Volume 47 Number 4 August 2006

Ferlat et al O-02.indd 441

441

18/12/2006 17:20:59

Proc. Fifth int. conf. on borate glasses, crystals and melts

application to water). These considerations motivated us to undertake an extensive description of B2O3 disordered phases with the SIESTA code. The purpose of the present work is to assess the validity of this approach by undertaking calculations on a series of known crystalline borate systems. To this end, results obtained with different numerical atomic orbitals (NAO) basis sets are confronted with calculations obtained using a plane-wave (PW) basis as well with experimental data. This allowed us to obtain efficient NAO basis sets to be used in forthcoming ab-initio molecular dynamics simulations of liquid and glassy B2O3.

2. Calculations 2.1. Technical details

First principles calculations were performed within the density functional theory (DFT) framework using the generalised gradient approximation (GGA) PBE exchange correlation functional.(23) Core electrons are replaced by norm-conserving Troullier–Martins pseudopotentials,(24) in the Kleinman–Bylander(25) form. Core radii of 1·4 and 1·7 a.u. were used for oxygen and boron, respectively. For calculations on B(OH)3 and B2O3 crystals, a 2×2×2 Monkhorst–Pack points(26) grid was used to sample the Brillouin zone. PW calculations on B(OH)3 were carried out with the PWSCF code(27) using a 50 Ry energy cutoff for the plane-wave expansion of the electronic wavefunctions. For SIESTA(20) calculations, the real-space integrations were carried out on a grid whose fineness is defined by a 200 Ry energy cutoff. The basis sets used in SIESTA are made of localised atomic orbitals: we used DZP (double-ζ polarised) basis sets which represent a good compromise between efficiency and accuracy. They amount to two 2s orbitals, two 2p shells and one 3d shell for both oxygen and boron atoms.

2.2. Generation of the NAO basis sets The performance of several DZP basis sets were tested and we present here the results obtained for two of them. These two basis differ in the cutoff radii of the support regions of their wave functions. The first basis, hereafter referred as DZPauto is obtained using the now standard procedure among SIESTA users: all radii are defined by one single parameter, the energy shift, i.e. the energy penalty payed by the orbital when confined.(28) This parameter was set to 150 meV. The second basis, hereafter referred as DZPopt is obtained using the variational basis optimisation scheme proposed in Anglada et al:(29) using a downhill-simplex minimisation, parameters defining the basis are varied so as to minimise

the enthalpy (E+PV) of a given reference system. The fictitious pressure, used to compress the orbitals radii and thus to reduce the computational efforts, was set to 0·2 GPa. The parameters that were allowed to vary were the orbitals cutoff radii and fractional charges for each atom type. The latter parameters influence the shape of an orbital: in cations, orbitals tend to shrink while they tend to swell for anions. Therefore, introducing a δQ parameter in the basis generating free atom calculations gives orbitals better adapted to ionic situations in the condensed system.(30) The reference system is a molecule or solid in which the atoms considered have a prominent role and which is small enough to allow many self consistent calculations with different basis parameters. We choose the boroxol molecule with hydrogens to cap the nonbridging oxygens as reference system (B3O6H3). The hydrogens pseudopotentials and basis parameters were taken from a recent study of liquid water.(22,31) The parameters used to define each basis set are summarised in Table 1. It is observed that the range of extension of the orbitals radii is about the same for both basis (in between 5 and 6 a.u. for the occupied ones) but while DZPauto provides cutoff radii longer for boron than for oxygen, DZPopt does the opposite. As will be shown in the next sections, the latter basis provides better structural results. In addition, the computational cost using DZPopt is reduced by about 50% as compared to calculations carried out with DZPauto. As for the boroxol ring molecule, both basis allowed to obtain interatomic distances which match the experimental values (obtained in the glass) and Hartree–Fock values(12) (obtained with a STO-3G basis) within 1%. The highest discrepancy (1·2%) is found for the B–O intra-ring distance (1·386 versus 1·37 Å). We shall now test the transferability and the performances of the obtained DZP basis by carrying out calculations on boric acid and boron trioxide crystals.

2.3. Boric acid B(OH)3 Boric acid B(OH)3, also known as sassolite mineral, is a molecular crystal having a layered structure with hydrogen bonding within the layers. It was determined by x-ray diffraction(32) that the crystals are – triclinic with four molecules in the unit cell and P1 space group symmetry. Periodic Hartree–Fock Table 1. Parameters describing the DZP basis sets rcB (a. u.) p d

rcO (a. u.) p d

δQB |e|

δQO |e|

s

DZPauto DZPopt

5·16 6·30 6·30 3·56 4·24 4·24 0 0 4·52 5·54 4·33 5·43 6·30 3·84 1·60 −0·16

s

442 Physics and Chemistry of Glasses: European Journal of Glass Science and Technology Part B Volume 47 Number 4 August 2006


18/12/2006 17:20:59


calculations have been carried out previously.(33) Using the experimental structure,(32) we performed relaxations of the coordinates using a PW basis and the two DZP basis derived in the previous subsection. The results are compared in Table 2. We note that DZP results (for either basis) compare very well with the PW ones. This means that the errors due to the use of localised basis sets is very small. It is indeed smaller that the errors introduced by the use of pseudopotentials and/or the use of the PBE exchange-correlation functional (as seen from the comparison of the PW results with the experimental data). Errors arising from whichever calculation on the B–O distance are about 2% (errors involving hydrogen atoms are much bigger due to the use of a pseudopotential for this atom but since our aim is the study of B2O3 we are not concerned by these errors).

2.4. B2O3 crystals We now turn to the description of B2O3 crystals. There are two known polymorphs, hereafter referred as B2O3-I and B2O3-II, in which the boron atoms have different coordination numbers(34,35) (see for instance Figures 1 and 2 in Takada et al(15)). B2O3-I is stable at ambient pressure while B2O3-II is a high pressure phase. The former has hexagonal symmetries while the latter is orthorhombic. A former ab-initio (plane wave pseudopotential) study carried out within the local density approximation (LDA) did not succeed in reproducing the experimental energy ordering and phase transition pressure between the two phases.(15) It was later shown that this failure was due to the use of the LDA: a proper account could be obtained(17) using the Perdew & Wang’s GGA.(36) We shall use the more recent PBE GGA functionality and show that our localised basis sets can well reproduce the transition features as well. The equilibrium curves for each crystal were obtained by allowing the atomic coordinates to relax from the experimental ones for different values of the volume cell. As in previous studies,(15,17) we considered only isotropic volume changes. The data obtained with the DZPopt basis are shown Figure 1 (very similar curves were obtained with DZPauto). The values for the equilibrium volume (relative Table 2. Interatomic distances (in Å) obtained in B(OH)3 crystal using different basis sets

DZPauto DZPopt PW (50 Ry)

Expt, Ref. 32

B–O O1–H1 O2–H2 O3–H3 B–Bintralayer B–Binterlayer

1·381 1·01 1·01 1·01 4·08 7·64

1·36 0·96 0·83 0·91 4·09 7·66

1·376 1·00 1·00 1·00 4·08 7·65

1·387 1·01 1·01 1·01 4·08 7·64

Figure 1. Energy (per B2O3 unit) as a function of the volume cell for B2O3-I and B2O3-II, crystals. The continuous lines are polynomial fits to the data and provided as guide for the eye. The values have been shifted to make the minimum of the B2O3-I curve equal to zero. The dashed straight line indicates the tangent line to both equilibrium curves (Gibbs construction): its slope is used to determine the transition pressure (4·1 GPa) to the experimental ones) are reported in Table 3 as well as the energy and volume difference between the two polymorphs. A slight expansion of the equilibrium volume cell is obtained (of the order of 5 and 2% for B2O3-I and B2O3-II, respectively) which results from the tendency of GGA functionals to overestimate the interatomic bond lengths. However, the B2O3-I crystal appears more stable than the B2O3-II, one, in agreement with experimental observations. Moreover, the values for the pressure transition between the two polymorphs (obtained from the Gibbs construction, see Figure 1) agree fairly well with the somehow uncertain experimental determination of 2 GPa.(37,38) We note however a dependency on the obtained value with the basis used (see Table 3). It is worthwhile to note that by reducing the B2O3-I volume cell to 73% of its original value, a transition to a new structure is observed in which part of the boron atoms became fourfold coordinated (and part of the oxygens became threefold coordinated). This Table 3. Equilibrium volumes (Veq) and energies (Eeq) obtained in B2O3 crystals using different basis sets. Vexp stand for the experimental known equilibrium volumes.(34,35) The transition pressures are obtained from the Gibbs construction (slope of the tangential line to the equilibrium curves)

DZPauto

DZPopt

PW (Engberg(17))

B2O3-I: (Veq−Vexp)/Vexp (%) 4·7 5·5 5·4 B2O3-II: (Veq−Vexp)/Vexp (%) 1·5 2·7 2·0 Veq(B2O3-II)-Veq(B2O3-I) (bohr) 61 65 54·8 Eeq(B2O3-II)-Eeq(B2O3-I) (mHa) 2·65 7·8 5·8 Ptransition (GPa) 1·4 4·1 3·1

Physics and Chemistry of Glasses: European Journal of Glass Science and Technology Part B Volume 47 Number 4 August 2006


443

18/12/2006 17:21:00


is a pressure induced transition although the final structure is not exactly the same as B2O3-II. Both phases are found to be insulators: electronic gap of 7·1 and 9·1 eV are obtained for B2O3-I and B2O3-II: respectively (a previous study, (16) carried out within the LDA, gave values of 6·2 and 8.85 eV, respectively).

3. Conclusions Ab-initio calculations using the SIESTA code have been carried out on a series of borate systems: the boroxol ring molecule: the boric acid crystal and the two known B2O3 crystals. Special efforts were devoted to obtaining and testing different atomic orbitals basis sets by comparing the obtained results to plane-wave calculations and experimental data. In all cases: it is shown that results of comparable accuracy to PW ones can be obtained with a rather small basis (DZP like). Indeed: the errors introduced by the use of localised orbitals appeared to be smaller than other sources of errors in the calculations (such as the use of the PBE exchange-correlation functional and pseudopotentials). The magnitude of these errors were found to be within the typical error-bar of DFTGGA calculations (i. e. ~2% on interatomic distances: ~5% on volume cell). The derived basis sets have been used for molecular dynamics simulations of B2O3 liquid and the results will be reported elsewhere. Our aim is an extensive exploration of the liquid at high pressure and high temperature as well as the description of the amorphous phases. We note that the characteristics of the high pressure transition between B2O3 polymorph are well reproduced by our calculations: a value of about 4 GPa is obtained for the transition pressure and the change of boron from three- to fourfold coordination can be reproduced by the compression of B2O3-I. This is particularely encouraging for the on-going high pressure studies of the liquid and amorphous phases.

Acknowledgements Technical support from A. P. Seitsonen is acknowledged. Part of the calculations have been performed

at the French supercomputing center IDRIS. G.F. would like to thank M. Fernandez-Serra for interesting discussions.

References

1. Nicholas, J., Sinogeikin, S., Kieffer, J. & Bass, J. Phys. Rev. Lett., 2004, 92, 215701. 2. Lee, S. K., Mibe, K., Fei, Y., Cody, G. D. & Mysen, B. O. Phys. Rev. Lett., 2005, 94, 165507. 3. Brazhkin, V., Katayama, Y., Inamura, Y., Kondrin, M., Lyapin, A., Popova, S. & Voloshin, R., JETP Lett., 2003, 78, 393. See also erratum in JETP Lett., 2004, 79, 308. 4. Hannon, A. C., Grimley, D. I., Hulme, R. A., Wright, A. C. & Sinclair, R. N. J. NonCryst. Solids, 1994, 177, 299. 5. Swenson, J. & Borjesson, L. Phys. Rev. B, 1997, 55, 11138. 6. Soules, T. F. J. Chem. Phys., 1979, 71, 4570. 7. J. Chem. Phys., 1980, 73, 4032. 8. Soules, T. F. & Varshneya, A. K., J. Am. Ceram. Soc., 1981, 64, 145. 9. Inoue, H., Aoki, N. & Yasui, I., J. Am. Ceram. Soc., 1987, 70, 622. 10. Takada, A., Catlow, C. R. A. & Price, G. D. J. Phys., Condens. Matter, 1995, 7, 8693. 11. Maranas, J. K., Chen, Y., Stillinger, D. K. & Stillinger, F. H. J. Chem. Phys., 2001, 115, 6578. 12. Zhuang, H.-Z., Zou, X.-W., Jin, Z.-Z. & Tian, D.-C. Phys. Rev. B, 1995, 52, 829. 13. Uchino, T. & Yoko, T. J. Chem. Phys., 1996, 105 (10), 4140. 14. Zwanziger, J. W. Solid State Nucl. Magn. Reson., 2005, 27, 5. 15. Takada, A., Catlow, C. R. A., Lin, J. S., Price, G. D., Lee, M. H., Milman, V. & Payne, M. C. Phys. Rev. B, 1995, 51, 1447. 16. Li, D. & Ching, W. Y. Phys. Rev. B, 1996, 54, 13616. 17. Engberg, U. Phys. Rev. B, 1997, 55, 2824. 18. Teter, M. Ab-initio based studies of borate glasses. Proc Second Int. Conf on Borates Glasses, Crystals and Melts. Soc. Glass Tech., Sheffield, 1997, 407–16. 19. Cormack, A. N. & Park, B. Phys. Chem. Glasses, 2000, 41 (5), 272–7. 20. Soler, J. M., Artacho, E., Gale, J. D., Garcia, A., Junquera, J., Ordejon, P. & Sanchez Portal, D. J. Phys., Condens. Matter, 2002, 14, 2745. 21. Kim, J., Mauri, F. & Galli, G. Phys. Rev. B, 1995, 52, 1640. 22. Fernández-Serra, M. V. & Artacho, E. J. Chem. Phys., 2004, 121, 11136. 23. Perdew, J. P., Burke, K. & Ernzerhof, M. Phys. Rev. Lett., 1996, 77, 3865. 24. Troullier, N. & Martins, J. L. Phys. Rev. B, 1991, 43, 1993. 25. Kleinman, L. & Bylander, D. M. Phys. Rev. Lett., 1982, 48, 1425. 26. Monkhorst, H. J. & Pack, J. D. Phys. Rev. B, 1976, 13, 5188. 27. Baroni, S. PWSCF. Http://www.pwscf.org. 28. Artacho, E., Sanchez-Portal, D., Ordejon, P. & Garcia, A., Phys. Status Solidi B, 1999, 215, 809. 29. Anglada, E., Soler, J. M., Junquera, J. & Artacho, E. Phys. Rev. B, 2002, 66, 205101. 30. Junquera, J., Paz, O., Sanchez-Portal, D. & Artacho, E. Phys. Rev. B, 2001, 64, 235111. 31. Fernández-Serra, M. V., Ferlat, G. & Artacho, E., Mol. Simul., 2005, 31, 361. 32. Zachariasen, W. H. Acta Crystallogr., 1954, 7, 305. 33. Zapol, P., Curtiss, L. A. & Erdemir, A. J. Chem. Phys., 2000, 113, 3338. 34. Gurr, G. E., Montgomery, P. W., Knutson, C. D. & Gorres, B. T. Acta Crystallogr., 1970, B26, 906. 35. Prewitt, C. T. & Shannon, R. D., Acta Crystallogr, 1968, B24, 869. 36. Perdew, J. P., Burke, K. & Wang, Y. JOURNAL TITLE? 1991. 37. Dachille, F. & Roy, R. J. Am. Ceram. Soc., 1959, 42, 78. 38. Mackenzie, J. D. & Claussen, W. F. J. Am. Ceram. Soc., 1961, 44, 79.

444 Physics and Chemistry of Glasses: European Journal of Glass Science and Technology Part B Volume 47 Number 4 August 2006


18/12/2006 17:21:00