Improved DFT calculation of Raman spectra of silicates

Vibrational Spectroscopy 56 (2011) 265–272

Contents lists available at ScienceDirect

Vibrational Spectroscopy journal homepage: www.elsevier.com/locate/vibspec

Improved DFT calculation of Raman spectra of silicates Daniel M. Többens ∗ , Volker Kahlenberg Institute of Mineralogy and Petrography, University of Innsbruck, Innsbruck, Austria

a r t i c l e

i n f o

Article history: Received 17 January 2011 Received in revised form 14 March 2011 Accepted 1 April 2011 Available online 9 April 2011 Keywords: DFT calculation Hybrid functional Vibrational frequency calculation Silicates

a b s t r a c t Very accurate vibrational spectra of silicates are obtained from DFT calculations if the appropriate Hamiltonian is used. Theoretical considerations suggest that the Hartree–Fock component of ACM1 hybrid functionals should be 1/6 instead of 1/4 for this class of compounds. When applied to the PBE functional this removes the scaling error of the calculated vibrational frequencies. Calculations using this PBE(n = 6) functional in combination with optimized Gaussian basis sets result in very small remaining deviations between observed and calculated Raman shifts, with standard uncertainties of ≈3.5 cm−1 , maximum deviations of ≈10 cm−1 , and no significant systematic trends. This has been confirmed for a wide range of silicate structures, for which high-quality Raman spectra have been published: forsterite ␣Mg2 SiO4 (nesosilicate), ␥-Y2 Si2 O7 (sorosilicate), K2 Ca3 Si3 O10 (oligosilicate), K2 Ca4 Si8 O21 (phyllosilicate), and ␣-quartz SiO2 (tectosilictae). © 2011 Elsevier B.V. All rights reserved.

1. Introduction With advances in modern technology, Raman spectrometers have become smaller, cheaper, and easier to use. As a result, Raman spectroscopy is nowadays a standard method in the characterization of materials. Beyond simple fingerprinting, vibrational analysis is extremely useful for obtaining information about structural features of molecules and solids [1]. The theoretical tools used for vibrational analysis have also advanced very much in recent years. Among the most powerful is certainly the ab initio calculation of vibrational spectra from the structure alone, using quantummechanical methods [2]. The calculation of vibrational spectra of molecular systems is a long-established procedure, implemented in most of the relevant computer codes [3,4]. Until recently, the situation was different for crystalline systems. Only a few ab initio codes did permit the calculation of vibrational spectra of crystalline compounds. This has changed with the advent of comparatively user friendly software [5], combined with the spectacular increase in computational power giving rise to a growing number of publications dealing with structures of increasing complexity. Many of these deal with compounds from the mineralogically and technically important class of the silicates. The benefit of this method is not so much in the calculated frequencies themselves, which after all are easily gained from the experimental spectrum. But the calculations yield additional information about the mode symmetries and the movements of the atoms involved in the particular vibrations. These can be used to

∗ Corresponding author. Tel.: +43 512 507 5532. E-mail address: [email protected] (D.M. Többens). 0924-2031/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.vibspec.2011.04.002

understand changes in frequencies or other features of the experimental spectrum, like peak intensity or broadening that result from changes in pressure, temperature, or chemical composition. Furthermore, deviations from the ideal crystal structure and its expected spectrum can be identified, and the presence of impurity phases can be spotted. In order to join the information from the experimental spectrum with the information from the calculation, it is necessary to assign the observed bands to specific calculated modes. The foremost criterion used in this context is the agreement of the observed and calculated frequencies. To allow for an unequivocal assignment, the uncertainty of both values needs to be lower than the distance between neighboring modes. For moderately complex silicates, as presented herein, this is ≈20 cm−1 on average, and considerably less in the some regions of the spectrum. Experimentally, this level of accuracy does not present a problem. Modern Raman spectrometers reach spectral resolutions of 2 cm−1 and accuracies of 0.5 cm−1 [6]. On the theoretical side the accuracy is lower, but the necessary quality is approached. A recent publication [7] reported average deviations of 6 cm−1 and 7 cm−1 for forsterite and ␣-quartz, respectively, with maximum deviations of 22 cm−1 and 18 cm−1 , using models of practical complexity and the well-established B3LYP hybrid functional. Our group found similar accuracies, 8 cm−1 and 20 cm−1 , in calculations using the less sophisticated LDA/VWN local density approximation [8]. The more modern PBE0 hybrid functional [9] is known to yield very accurate geometries, but overestimates vibrational frequencies [7]. The agreement with experimental values is unimpressive. However, when we tested this functional on one of our silicate structures [10] we found that the deviations were highly systematic, with a direct proportionality

266

D.M. Többens, V. Kahlenberg / Vibrational Spectroscopy 56 (2011) 265–272

between calculated and experimental frequencies. After correcting the calculated values with an ad hoc scaling factor, the agreement was extremely good. This result directly triggered the formulation of this paper, in which we establish that the scaling factors for frequency calculations using PBE0 or the related PBE functional [11] are the same for all silicates. Furthermore, while the scaling of theoretical frequencies is an established approach [12–14], it is also unsophisticated, hinders subsequent thermodynamic calculations of enthalpy and entropy, and reduces the confidence that other researches will put in the results. However, the scaling factor was found to be a function of the fraction of Hartree–Fock exchange energy used in the hybrid functional. By selecting a Hartree–Fock component of 1/6, instead of the 1/4 used in the PBE0 functional, the scaling error is removed, and the calculated values can be used unmodified. This is in agreement with the underlying theory of the ACM1 functionals, which predicts the value of 1/6 for silicates and similar compounds.

2. Computational parameters All calculations were done from first principles using 3Dperiodic density functional theory with Gaussian basis sets implemented in the program CRYSTAL09 [15]. Default settings of the program were used unless otherwise mentioned. The vibrational frequencies were calculated in harmonic approximation at the -point [5], at a stationary point on the potential energy surface. Therefore, the structure was fully relaxed by symmetry-preserving but otherwise unconstrained geometry optimization of atomic positions and lattice parameters. First and second derivatives of the energy for the calculation of the mass weighted Hessian matrix were computed from small displacements of the atoms. The default step size of 0.001 au was used for the calculation of the numerical gradient; significant changes >1 cm−1 of the calculated frequencies were observed only for step sizes >0.005 au. A Pack–Monkhorst k net with 4 × 4 × 4 points in the Brillouin zone was employed for all calculations with the exception of forsterite and quartz, for which 6 × 6 × 6 points were used. With these values, the energy gain from using denser nets was <10−6 Eh and <10−7 %, which we found to be the limit of improvement for the accuracy of calculated vibrational frequencies. The level of numerical accuracy was increased over the default settings of the software as described i.e. in [16], selecting tolerances for coulomb and exchange sums with keyword TOLINTEG 7 × 7 × 7 × 7 × 15, thresholds of the SCF energy of 10−9 Eh for the geometry optimizations and of 10−10 Eh for the frequency calculations, as well as a (75,974)p grid (keyword XLGRID) for the numerical integration of the DFT exchange-correlation contribution. Modified Broyden damping was applied [17]. The basis sets of the atoms were selected as follows: for silicon a 86-311G(1) contraction [18], for oxygen a 8-411G(1) contraction based on the 8-411G contraction given in [19], for calcium a 86-511G(21) contraction used for grossular in [20], for potassium a 86-511G(21) contraction based on the 86-511G(3) contraction used in [21], for magnesium a 8-511G(1) contraction based on the 8-511G contraction given in [22]. Only for yttrium no allelectron calculation was conducted, but a Hay and Wadt small core – effective core pseudopotential was used in a HAYWSC-311G(11) contraction based on the unpublished expansion of Gennard and Furio given in the CRYSTAL basis sets library [23]. For all basis sets, the exponents of the two most diffuse sp shells and of the one most diffuse d shell have been re-optimized in the respective structures. The values are available as supporting information together with the optimized structures and frequency calculation files. For the vibrational spectra of quartz, all of the E and A2 modes are split at the zone center into transverse and longitudinal components. Since ␣-quartz is non-centrosymmetric, this TO-LO splitting

is also observed in the Raman spectrum, for the Raman-active E modes. The values of these splittings were calculated with the Berry phase approach, using the dielectric tensor ε11 = ε22 = 2.356, ε33 = 2.383 [24]. 3. Theory The PBE0 functional [9] is a hybrid functional based on the gradient corrected Perdew–Burke–Ernzerhof PBE functional [11]. It is known to give very exact geometries, which accounts for its popularity, but was found to perform unsatisfactory in the calculation of vibrational spectra of silicates [7]. It is shown in here that this is due to constraints used in the formulation of PBE0, which are correct for light-atom molecules, but not for solids based on heavy atoms. By choosing the correct conditions the PBE0 functional can be modified easily to give very accurate vibrational spectra for these compounds. The basic idea of hybrid functionals is to improve the density functional approximation for the exchange-correlation energy (DF) by replacing an appropriate fraction of the exchange energy ExDF with Hartree–Fock exchange energy ExHF . Within the family of functionals known as ACM1 the HF/DF exchange ratio is ruled by only one system- and property-dependent coefficient: hybrid

Exc

DF = Exc + a1 (ExHF − ExDF )

Perdew et al. demonstrated [25] that the optimum value of the a1 coefficient can be fixed a priori: hybrid

Exc

DF = Exc +

1 HF (E − ExDC ) n x

In this equation, the integer n describes the shape of energy as a function of the coupling constant. The optimum value of n is given by the lowest order of perturbation theory that provides a realistic description of this shape, and in most cases can be estimated by examining the convergence of the Møller–Plesset perturbation expansion. Since fourth-order perturbation theory (MP4) was found to be sufficient to get accurate numerical results for most molecular systems [26], n = 4 was suggested as the best single choice. The resulting family of hybrid functionals is known as ACM0. In case that for the DFT part of the calculations the generalized gradient approximation functional PBE is used, the PBE0 functional results: PBE0 PBE Exc = Exc +

1 HF (E − ExPBE ) 4 x

However, it was shown [27] that the convergence behavior of the MPn series is different for systems for which electron pairs (core, bond, lone pairs) are well separated and equally distributed over the whole space of the system in question (class A), and systems that are characterized by electron clustering in a confined region of space (class B). For the class A systems, the convergence is monotonic and fourth-order calculations are sufficient. Since this class contains (light) electropositive atoms and saturated hydrocarbons, n = 4 is indeed the correct choice for most systems composed of organic molecules. On the other hand, for class B systems, oscillations between even and odd orders of MP theory require calculations up to MP6 to obtain a reliable estimate of the correlation energy at infinite order [28]. Consequently for these systems the functional should be reformulated with n = 6, so that PBE(n=6)

Exc

PBE = Exc +

1 HF (E − ExPBE ) 6 x

Silicates, which are dominated by the electronegative oxygen, can be expected to be class B systems [28]. They should therefore be calculated with n = 6, hence with a modified functional called PBE

D.M. Többens, V. Kahlenberg / Vibrational Spectroscopy 56 (2011) 265–272

267

(n = 6) in the following. It should be noted that a similar hybridization factor of 0.16 has also been introduced recently, on an empirical basis, in the WC1LYP [7] and B1WC [29] Hamiltonians. Vibrational frequencies are calculated from a mass weighted Hessian matrix constructed from first and second derivatives of the energy differences that result from small displacements of the atoms. This procedure is obviously very sensitive to the shape of the energy function. 4. Performance It is demonstrated in the following that the approach with n = 6 indeed removes systematic errors from the calculated vibrational spectra of all silicates. 4.1. Forsterite Forsterite (␣-Mg2 SiO4 ) is an inosilicate. The structure consists of isolated SiO4 tetrahedra within a framework of edgesharing MgO6 octahedra. Forsterite has Z = 4 formula units in the unit cell, space group Pbnm (point group mmm or D2h ), and symmetry representation analysis gives the species of the optical modes as  = 11Ag (R) + 10Au (S) + 11B1g (R) + 9B1u (I) + 7B2g (R) + 13B2u (I) + 7B3g (R) + 13B3u (I). R, I, and S denote Raman and infrared activity and silent modes, respectively. Because of its small unit cell (≈290 A˚ 3 ) and low number of six symmetrically independent atoms DFT calculations of its properties are cheap. It is thus popular as an archetype of silicates in general. In recent times, it has been used as a test case for the DFT-based calculation of vibrational spectra at least twice: A calculation with the B3LYP functional [16] compared the results to four independent experimental Raman and IR spectra taken from literature [30–33]. Another work [7] used experimental Raman [30] and IR [34] spectra taken from the literature, together with other silicates, to gauge the performance of various established functionals in the calculation of vibrational spectra. In here, the experimental Raman shifts of Chopleas [30] were used (Table 1). These polarized single crystal Raman spectra come with experimentally determined mode symmetries for all 36 Raman active modes. Matching this information with the symmetries of the calculated frequencies minimizes the risk of cross-assignment. The uncertainties of the frequencies are given as ±1 cm−1 . The experimental conditions seem most suited to produce reliable data. As noted in [16], for some modes the frequencies given in the literature vary significantly (very broad or weak modes, or ones close to much stronger lines). In here, only those modes were considered reliable, for which the range of experimental values given by the four authors [30–33] is <5 cm−1 , and which have been observed by at least three of these authors (Table 1). Only these reliable values were used to calculate the statistical agreement indicators. In Table 1 the experimental data are compared with calculated Raman modes from various established functionals. The widespread B3LYP functional has been applied successfully to a wide variety of problems. In comparative calculations of vibrational spectra it was repeatedly found to give the best results [7,35]. It can be considered the current benchmark. The local density approximation, combining the Dirac–Slater (LDA) exchange functional [36] with the Vosko–Wilk–Nusair (VWN) correlation functional [37], has the advantage of being very fast, about four times faster than B3LYP. It was applied successfully to the Raman spectra of various silicates by our group. In some cases the results were superior to those derived from B3LYP [8]. The WC1LYP functional has been proposed recently [7] in the context of calculating vibrational spectra of silicates, and rivals B3LYP in quality of the results. It uses a hybridization factor of 0.16. The PBE0 functional is known to give

Fig. 1. Deviations ((obs) − (calc)) between observed [30] and calculated Raman frequencies of forsterite. Filled symbols designate reliable observed frequencies, open symbols unreliable values.

very exact geometries, but was found to perform not as good as B3LYP in the calculation of vibrational spectra [7]. PBE, the gradientcorrected functional the PBE0 hybrid functional is based on, was also calculated. By itself this functional does not perform well in the calculation of vibrations [7,35]. From the simple statistical indicators in Table 1 the calculations just re-affirm what has been stated before: WC1LYP gives the best results, followed by B3LYP and PBE0 performing with comparable quality, LDA/VWN is worse, and PBE-calculated frequencies are very bad. However, there is a significant difference between the deviations resulting from PBE-based Raman spectra (PBE, PBE0) and the others. In the latter cases the positive and negative deviations are scattered over the range of observed vibrational frequencies. In the case of B3LYP the strongest absolute deviations are located in the center of the spectrum, at ≈600 cm−1 . For LDA/VWN they are located at the boundaries of the spectrum, at ≈300 cm−1 and ≈1000 cm−1 . This cannot be corrected by a simple adjustment. The situation is entirely different for PBE and PBE0 (Fig. 1). It is immediately obvious that the major part of the deviation is systematic. Furthermore, the unsystematic part of the deviations that would remain after linear scaling is very similar for the two functionals. If the correct exchange ratio of 1/6 is used with n = 6, the systemic deviation vanishes. The results from this calculation, designated as PBE(n = 6) in the following, are given in Table 1 and Fig. 1. The remaining deviations from the observed values give agreement indicators of  = 2.2 cm−1 and ||max = 4.8 cm−1 , far better than any of the other functionals. The results indicate that the ideal value for the HF fraction would be slightly lower than the 1/6 value prescribed by theory, at a1 = 0.156. However, this small difference is not significant, resulting in changes <3 cm−1 for high-frequency modes. As it turned out this is true in all tested cases, with ideal values of a1 varying only slightly in the range 0.15–0.17. Besides the original PBE functional [11] a number of revisions have been published. Among these the one known as PBEsol [38] has been geared specifically towards improving “equilibrium properties of densely packed solids and their surfaces”, and has become quite popular. When this revised functional is used, the systematic trend of the deviations is reduced in comparison to the original PBE functional. Insofar the claim of improved description of solidstate properties is valid for the calculation of vibrations, too, even though these are not strictly equilibrium properties. But the results are inferior as far as the unsystematic deviations are concerned. This situation cannot be improved by scaling or hybridization. For this reason the use of the unmodified PBE functional is preferable for the calculation of vibrational spectra.

268

D.M. Többens, V. Kahlenberg / Vibrational Spectroscopy 56 (2011) 265–272

Table 1 Observed [30] and calculated Raman shifts [cm−1 ] for forsterite, matching modes with identical symmetry. Statistical indices of the differences between observed and calculated values are calculated from reliable values only (see Section 4.1). Observed

Reliable

175 183 183 – 220 226 226 242 242 274 286 304 304 315 315 318 318 323 323 329 329 339 339 351 365 365 374 374 383 410 422 422 434 434 439 545 545 582 586 586 592 608 608 632 632 824 824 838 838 856 856 866 866 881 881 920 920 965 965 975 975 Statistical indices, cm−1 Standard deviation  Largest negative deviation min Largest positive deviation max

Sym

B3LYP

LDA/VWN

WC1LYP

PBE

PBE0

PBE(n = 6)

PBEsol

B2g A1g B3g B1g A1g B2g B1g B3g A1g B3g B1g B2g A1g A1g B1g B2g B3g B1g B3g A1g B1g B2g A1g B1g B2g B3g A1g B1g A1g B1g A1g B1g B2g B3g A1g B1g

181.6 185.9 184.1 222.9 231.2 241.9 255.9 293.7 302.9 320.4 313.2 319.6 324.3 343.0 360.1 371.0 375.2 388.9 416.8 418.4 436.0 442.7 559.9 597.9 603.2 605.9 617.8 642.3 826.0 841.7 861.5 871.5 887.3 930.2 968.6 981.0

178.5 192.3 200.3 228.3 227.0 252.3 263.2 288.7 325.4 331.0 331.9 338.1 339.7 364.5 369.9 377.9 391.5 391.9 422.4 440.5 451.0 453.3 539.1 583.4 579.8 588.9 608.2 632.9 830.4 843.9 859.2 869.4 889.6 928.2 981.6 992.1

179.7 186.6 184.7 223.5 229.1 242.4 256.5 289.6 308.7 317.5 322.3 326.4 330.2 345.9 359.7 370.6 380.9 388.4 414.2 422.1 439.0 441.1 551.5 590.4 595.5 600.7 612.8 637.1 825.7 840.1 859.0 868.6 888.9 930.1 971.8 983.7

173.0 181.0 181.9 213.7 220.2 235.3 247.7 279.5 297.1 309.1 306.1 311.1 316.5 331.5 346.6 357.7 363.7 376.5 400.0 407.8 423.5 425.5 525.5 561.2 566.7 571.9 589.0 610.7 793.6 807.0 821.4 831.5 844.7 885.4 934.8 945.4

180.3 189.7 194.4 224.8 229.8 246.9 259.3 290.4 313.6 323.0 323.9 325.2 333.6 348.5 364.6 373.7 380.4 392.4 419.3 427.6 442.6 446.4 556.3 596.2 599.1 604.2 619.4 643.7 843.4 857.4 870.7 881.1 897.3 938.8 983.1 994.5

177.9 186.6 189.9 221.0 226.6 242.6 254.4 286.5 307.4 317.4 318.5 319.7 327.3 342.7 358.9 368.8 374.2 386.9 413.4 420.5 435.4 439.3 546.6 584.8 588.7 593.8 609.5 632.8 828.3 842.1 855.8 866.0 881.6 922.9 968.6 979.8

173.3 185.4 191.3 218.6 221.0 241.2 253.1 279.2 309.8 316.6 317.9 321.1 326.5 344.2 354.3 363.8 373.8 382.3 406.4 421.3 433.7 434.3 526.7 566.9 566.3 574.0 593.4 616.1 809.0 822.0 837.1 847.3 863.3 902.4 954.5 965.0

5.6 −17.2 4.8

8.2 −25.5 6.2

3.0 −10.1 0.5

11.3 2.0 36.3

5.5 −19.5 −2.2

2.2 −4.8 3.3

8.8 −5.8 19.7

4.2. Quartz Quartz (␣-SiO2 ) is at the other extreme of the range of silicate structure building principles. All SiO4 tetrahedra are cornerconnected in a three-dimensional framework. It is the prime example of the tectosilicates, though it should be mentioned that formally quartz is denoted as a tectosilicate only in the Dana classification system (75.01.03.01), but as an oxide in the Strunz system (04.DA.05). Quartz has Z = 3 formula units in the unit cell, space group P31 21 (point group 32 or D3 ), and symmetry representation analysis gives the species of the optical modes as  = 4A1 (R) + 4A2 (I) + 8E(R,I). All of the E modes are split into transverse and longitudinal components. Its qualities as a test case are similar to forsterite: (i) the unit cell is small and thus calculation costs are low, and (ii) multiple high-quality Raman spectra with experimentally determined symmetry modes are available in the literature [39–41]. Unsurprisingly, this system has been used to assert the quality of calculated vibrational spectra for different functionals and basis sets [7,35]. In the present study, the observed line positions from [39] were used (Table 2). The experimental uncertainties of these values are <1 cm−1 . In this work the E–E(LO) pair at 452/512 cm−1 has not been observed, as it is hidden by a strong A1 mode. These values were taken from [40] instead. Again, for some modes the frequencies given in the literature vary. In here, only those modes are considered reliable, which have been observed in all three papers [39–41] with frequencies within a range of 6 cm−1 .

Fig. 2. Deviations ((obs) − (calc)) between observed [39] and calculated Raman frequencies of ␣-quartz. Filled symbols designate reliable observed frequencies, open symbols unreliable values.

Table 2 compares the calculated mode frequencies with the observed values. The situation is the same as observed for forsterite (Fig. 2): The larger part of the deviations in the calculations using the PBE and PBE0 Hamiltonians are systematic, and these systematic deviations are of the same magnitude. The agreement of the PBE(n = 6) calculation with the observed frequencies is slightly worse than in the case of forsterite, with

D.M. Többens, V. Kahlenberg / Vibrational Spectroscopy 56 (2011) 265–272

269

Table 2 Observed [39] and calculated Raman shifts [cm−1 ] for ␣-quartz, matching modes with identical symmetry. Statistical indices of the differences between observed and calculated values are calculated from reliable values only (see Section 4.2). The PBE(n = 6)d(11) frequencies have been calculated with the extended basis set described in Section 4.2. Observed

Reliable

128.0 128.0 205.6 263.1 263.1 354.3 354.3 393.8 393.8 452.5 [40] 463.6 463.6 697.4 697.4 796.7 796.7 1066.1 1066.1 1083.0 1083.0 1160.6 1160.6 128.0 128.0 263.1 263.1 401.8 401.8 512.0 [40] 697.4 697.4 808.6 808.6 1231.9 1160.6 Statistical indices, cm−1 Standard deviation  Largest negative deviation min Largest positive deviation max

Sym

B3LYP

WC1LYP

PBE

PBE0

PBE (n = 6)

PBE (n = 6) d(11)

PBEsol

E A1 E A1 E E A1 E E E A1 E E(LO) E(LO) E(LO) E(LO) E(LO) E(LO) E(LO) E(LO)

131.9 195.1 268.5 356.1 400.6 450.2 467.8 693.6 792.4 1074.5 1090.4 1173.4 131.9 270.3 408.8 512.6 695.4 806.6 1236.4 1171.2

133.3 212.0 267.6 350.6 395.7 450.1 469.5 697.8 798.0 1078.3 1093.8 1169.5 133.3 269.8 405.2 511.3 700.2 811.8 1241.3 1166.9

127.9 211.0 255.9 334.7 380.8 428.2 447.1 674.2 773.1 1036.3 1053.7 1127.7 127.9 258.2 388.4 489.0 676.5 785.0 1205.4 1125.5

135.2 218.5 269.4 353.4 399.7 454.7 473.6 704.4 805.8 1087.5 1103.7 1180.1 135.2 271.5 408.4 516.0 706.9 819.5 1250.3 1177.5

132.9 216.1 265.3 347.7 393.9 446.6 465.9 695.0 795.6 1071.8 1088.2 1164.0 132.9 267.5 402.4 507.7 697.5 808.7 1236.4 1161.6

132.4 208.7 264.4 347.1 392.6 442.0 462.0 693.3 794.8 1064.0 1080.9 1161.4 132.4 266.4 401.1 499.9 695.7 807.2 1230.4 1159.2

131.2 231.3 258.2 331.1 378.1 432.6 459.2 688.6 790.0 1051.9 1069.6 1134.0 131.2 261.3 387.8 491.7 692.0 801.7 1221.5 1131.5

4.9 −12.8 4.3

4.1 −12.2 3.7

10.5 0.1 32.9

6.1 −21.4 0.9

3.4 −5.7 6.6

3.1 −4.4 7.2

8.6 −3.2 26.6

For Raman modes not observed at all in [39] the positions given in [40] are listed. Longitudinal-optical modes (LO) are listed in the same order as the corresponding transversal-optical modes.

 = 3.4 cm−1 and ||max = 6.6 cm−1 . Simultaneously, the B3LYP Hamiltonian performs better than in the forsterite case, while the result from WC1LYP is not as good. However, the use of PBE(n = 6) still results in considerable improvement over the best available alternatives. A larger basis set was tested, using two d shells both for oxygen and silicon. This resulted in further increase of the quality, but not decisively so. This is in agreement with earlier studies on the influence of the basis set size on the calculated vibrational frequencies [35]. Notably, the longitudinal-optical modes are calculated with the same level of accuracy as the transversal ones. 5. Application to other silicates

functional (Fig. 3). The numerical values are available as supporting information. This new calculation confirmed the results of the previous work. All the mode assignments stayed the same as derived from the LDA/VWN model. The accuracy of the calculated frequencies increased over the older calculations, with agreement factors of  = 3.1 cm−1 , min = −10.4 cm−1 , max = 5.1 cm−1 . With this increased reliability of the mode assignments it becomes definite that various features of the observed Raman spectrum cannot be explained within the framework of the model calculations, that is they are not proper 1st order Raman modes of the 1st Brillouin zone of the ideal crystal. All of these features are located in the s(Si–O) range of the spectrum – two weak unexplained lines at 891 cm−1 and 999 cm−1 , and the apparent splitting of the lines at 1048 cm−1 and 1105 cm−1 into pairs. The most likely explanation

In the following, the PBE(n = 6) functional is applied to a number of other silicates, for which high-quality Raman spectra have been published. This allows to gauge the performance in praxis, with more complex structures. Additionally, different silicate groups with varying degrees of polymerization are tested. As is often the case in praxis, no experimental symmetry species are available. In the absence of computationally expensive intensity calculations, the accuracy of the calculated frequencies is crucial in avoiding cross-assignments. 5.1. K2 Ca4 Si8 O21 K2 Ca4 Si8 O21 is a single layer silicate [6]. It is triclinic, space group P1 (point group 1 or Ci ), with a unit cell volume of 481 A˚ 3 containing Z = 1 formula unit. The structure is unusual in that it contains both secondary (Q2 ) and tertiary (Q3 ) tetrahedra. Symmetry representation analysis gives the species of the optical modes as  = 51Ag (R) + 51Au (I). Experimentally, 44 Raman bands were observed. Originally, the Raman spectrum was interpreted using a pure DFT-calculation based on the LDA/VWN Hamiltonian. This approach already worked fairly well, with deviations of  = 5.9 cm−1 , min = −15.5 cm−1 , max = 8.0 cm−1 [6]. In here, the vibrational spectra have been recalculated using the PBE(n = 6)

Fig. 3. Experimental Raman spectrum of K2 Ca4 Si8 O21 , observed line positions (top row) and values calculated with LDA/VWN (bottom) taken from the literature [6]. The middle row of bars gives the line positions calculated with the PBE(n = 6) functional. The agreement is superior, but no significant changes in the assignment occur. Note that some weak lines and splitting in the observed spectrum remain unexplained.

270

D.M. Többens, V. Kahlenberg / Vibrational Spectroscopy 56 (2011) 265–272

Fig. 4. Experimental Raman spectrum of K2 Ca3 Si3 O10 , observed line positions (top row), and frequencies calculated with the PBE(n = 6) functional (middle row). The frequencies in the lower row were calculated using PBE0 (with n = 4) and required the application of an ad hoc scaling factor of 0.977 [10].

Fig. 5. Experimental Raman spectrum of ␥-Y2 Si2 O7 , observed line positions (top row) and values calculated with LDA/VWN (bottom). The middle row of hash marks gives the line positions calculated with the PBE(n = 6) functional. The increased quality is obvious.

is the presence of a small amount of an unknown 2nd phase; there exists a large number of different K–Ca silicates, and their phase regimes are still not fully understood. Even though the Raman spectrum was obtained from a single-crystal sample, it is possible that a small fraction from another phase is intergrown and escaped detection (Fig. 4).

tered unit cell, space group P21 /c (point group 2/m or C2h ), and symmetry representation analysis gives the species of the optical modes as  = 15Ag (R) + 17Au (S) + 15Bg (R) + 16Bu (I). The Raman spectrum (Fig. 5) was originally published and discussed as part of a comparative study [8], in which the local density approximation based on Dirac–Slater (LDA) exchange and Vosko–Wilk–Nusair (VWN) correlation functionals was employed. With this setting, the experimental Raman spectra of all polymorphs could be modelled with average deviations of 8 cm−1 and maximum deviations of ±20 cm−1 . While in the original comparative study yttrium was described by a HAYWSC-311G(1) contraction, in here we used a HAYWSC-311G(11) contraction that more closely conforms to the original unpublished expansion of Gennard and Furio given in the CRYSTAL basis sets library [23]. This resulted in some changes, but without decisive quality improvements. Of the total 30 Raman active modes some could be interpreted unequivocally, in particular the Si–O–Si bending mode at 661 cm−1 that is characteristic for the [Si2 O7 ] units. However, in the region of the Si–O stretching modes (900–1000 cm−1 ) the assignments remained arguable due to the close proximity. The same holds for the 400–450 cm−1 region, where Y–O stretching modes and O–Si–O bending modes concur. With the PBE(n = 6) calculation the agreement factors fall to  = 3.3 cm−1 , min = −4.7 cm−1 , max = 7.8 cm−1 . This allows the unequivocal assignment of previously doubtful modes, e.g. the mode at 433 cm−1 as O–Si–O bending.

5.2. K2 Ca3 Si3 O10 K2 Ca3 Si3 O10 [10] is an oligosilicate with trimeric basic building units in a strongly bent (horseshoe) conformation. Equally oriented [Si3 O9 ]-groups are arranged in a layer-like fashion, but remain isolated. K2 Ca3 Si3 O10 has Z = 2 formula units in the unit cell, space group C2/c (point group 2/m or C2h ), and symmetry representation analysis gives the species of the optical modes as  = 26Ag (R) + 25Au (S) + 28Bg (R) + 26Bu (I). Originally, in the work that prompted this paper, the interpretation of the observed Raman spectra was based on a PBE0 calculation. An ad hoc scaling factor of 0.977 was derived from unequivocal modes in the high-frequency range and was applied to the calculated frequencies. This resulted in a very good description of the observed spectrum, with agreement factors of  = 2.8 cm−1 , min = −5.9 cm−1 , max = 6.9 cm−1 , and unequivocal assignment of nearly all modes, both experimental and calculated, limited only by the experimental resolution and peak width. Compared to this the results from the PBE(n = 6) calculation are slightly inferior, with agreement factors of  = 3.4 cm−1 , min = −6.4 cm−1 , max = 8.5 cm−1 . However, no scaling was necessary, especially not one based on the analyzed compound itself. Given the high quality of both calculations, no significant differences in the interpretation are expected. In fact, a number of formal mode assignments change, where small changes in the calculated frequencies put a different mode closer to the observed band. However, all of these are insignificant. The experimental peaks are so wide as to cover the whole range of the modes involved, and most probably consist of multiple modes. 5.3. -Y2 Si2 O7 ␥-Y2 Si2 O7 [42] is a sorosilicate, composed of [Si2 O7 ] dimers. The two tetrahedra are in antiparallel (staggered) orientation. Due to the inversion symmetry of the [Si2 O7 ] moiety the Si–O–Si angle at the oxygen connecting the tetrahedra is 180◦ . The compound is a metastable polymorph, one of seven currently known forms with this composition. ␥-Y2 Si2 O7 has Z = 2 formula units in the cen-

6. Geometry While the focus of this article is on the calculation of vibrational spectra, other parameters cannot be ignored completely. Most notably the accuracy of the calculated crystal structure is of fundamental importance to any further analysis, including the vibrational analysis. Table 3 reports the relevant structural data for the PBE(n = 6) functional. The interatomic distances and lattice parameters are systematically overestimated by about 1%, the unit cell volume by 3%. These deviations are slightly higher than the overestimation resulting from the use of the PBE0 functional, which for the same structures and basis sets overestimates the unit cell volume by 2%. This is as expected, since the PBE functional strongly overestimates distances [7]. Lattice angles of the unit cell in the non-orthogonal compounds are excellently reproduced, with deviations less than 1◦ . The same is true for the Si–O–Si angles between the corner-connected tetrahedra, and the O–Si–O bond angles within the SiO4 tetrahedra are reproduced even more

D.M. Többens, V. Kahlenberg / Vibrational Spectroscopy 56 (2011) 265–272

271

Table 3 Agreement of calculated and experimental lattice parameters, interatomic distances and angles and cell volume when using the PBE(n = 6) functional (in Å or degree, respectively). Parameter Lattice parameters, Å Forsterite

␥-Y2 Si2 O7

K2 Ca3 Si3 O10

K2 Ca4 Si8 O21

␣-Quartz Lattice parameters, angles, deg ␥-Y2 Si2 O7 K2 Ca3 Si3 O10 K2 Ca4 Si8 O21

Unit cell volume, Å3 Forsterite ␥-Y2 Si2 O7 K2 Ca3 Si3 O10 K2 Ca4 Si8 O21 ␣-Quartz Average bond length, Å Forsterite ␥-Y2 Si2 O7 K2 Ca3 Si3 O10

K2 Ca4 Si8 O21

␣-Quartz Average bond angle, deg Forsterite ␥-Y2 Si2 O7 K2 Ca3 Si3 O10 K2 Ca4 Si8 O21 ␣-Quartz ␥-Y2 Si2 O7 K2 Ca3 Si3 O10 K2 Ca4 Si8 O21 ␣-Quartz

Exp. a b c a b c a b c a b c a c

Calc.

4.755 10.196 5.981 4.660 10.780 5.540 10.362 10.581 9.831 6.796 7.085 11.197 4.914 5.405

4.778 10.252 6.002 4.730 10.915 5.621 10.402 10.746 9.878 6.857 7.162 11.276 4.938 5.439

Dev. 0.023 0.056 0.021 0.070 0.135 0.081 0.040 0.165 0.046 0.061 0.077 0.079 0.024 0.034

ˇ ˇ ␣ ˇ 

96.10 118.23 96.69 105.34 109.23

95.98 117.51 96.63 105.38 109.29

−0.12 −0.72 −0.06 0.04 0.06

Volume Volume Volume Volume Volume

290.0 276.7 949.7 478.5 113.0

294.0 288.6 979.3 491.2 114.8

4.0 11.9 29.6 12.7 1.8

Mg[6]–O Si–O Y[6]–O Si–O K[10]–O Ca[6,8]–O Si–O K[8]–O Ca[6]–O Si–O Si–O O–Si–O O–Si–O O–Si–O O–Si–O O–Si–O Si–O–Si Si–O–Si Si–O–Si Si–O–Si

2.112 1.636 2.285 1.612 2.896 2.481 1.625 2.929 2.418 1.619 1.609 109.15 109.36 109.32 109.32 109.47 180.00 160.22 133.66 180.00 143.60

exact. The most notable structural deviation of all calculations was the O–Si–O angle in K2 Ca3 Si3 O10 , which was too low by −7◦ . This deviation persists (with −5◦ ) when the original PBE0 functional is used [10]. The high value in the real structure is experimentally confirmed, and is in agreement with systematic trends in trimeric silicates. The reason for this deviation is unknown. 7. Summary Summing up, the overall quality of the vibrational spectra calculated with the PBE(n = 6) functional is exceedingly good, with standard deviations of 2–4 cm−1 , maximum deviations of 10 cm−1 and below, and no significant scaling error. The HF/DF exchange ratio of 1/6 is the correct choice for silicates. The applicability of the PBE(n = 6) setup to compounds besides silicates is still an open question. Preliminary results, based on individual calculations, indicate that it also works very well for aluminates. The results

2.120 1.647 2.287 1.639 2.939 2.509 1.639 2.950 2.433 1.635 1.623 109.12 109.45 109.36 109.31 109.47 180.00 153.35 134.46 180.00 142.51

0.008 0.011 0.003 0.026 0.043 0.028 0.014 0.021 0.015 0.017 0.014 −0.03 0.09 0.04 −0.01 0.00 – −6.87 0.80 – −1.09

for corundum, ␣-Al2 O3 , are provided as supplementary information. The agreement of calculated frequencies of the seven Raman active modes with experimental data [43] is very good, with agreement factors of  = 2.6 cm−1 , min = −3.4 cm−1 , max = 4.5 cm−1 . For borates, our preliminary results indicate that the accuracy falls behind the B3LYP functional. This is in agreement with the theoretical background, as the difference in electronegativity is even larger for the Al–O combination than for Si–O, at otherwise very similar properties. Boron-dominated structures, on the other hand, can be expected to behave as class A systems [27]. If instead of modifying the hybridization factor scaling is used, the quality of the agreement does not suffer. For the frequencies calculated with the PBE functional, the best scaling factor, averaged from forsterite and ␣quartz, is 1.033. Using the pure DFT functional PBE is faster, which is useful especially for complicated structures. On the other hand, a further increase of the available computational power might allow incorporating the effects of anharmonicity and thermal expan-

272

D.M. Többens, V. Kahlenberg / Vibrational Spectroscopy 56 (2011) 265–272

sion into the calculations. While these are comparatively small for inorganic solids, with the current level of accuracy they are not negligible anymore. Finally, in this study only Raman spectra were used for their higher spectral resolution. Nevertheless, there is no reason why the calculated infrared frequencies should be of inferior quality. Acknowledgements This work was supported by the Austrian Ministry of Science BMWF as part of the UniInfrastrukturprogramm of the Forschungsplattform Scientific Computing at LFU Innsbruck. Financial support from Sandoz GmbH is gratefully acknowledged. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.vibspec.2011.04.002. References [1] J. Chalmers, P. Griffiths (Eds.), Handbook of Vibrational Spectroscopy, Vols. 4–5, J. Wiley & Sons, Chichester, 2001. [2] P. Pulay, G. Fogarasi, F. Pang, J.E. Boggs, J. Am. Chem. Soc. 101 (10) (1979) 2550–2560. [3] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery Jr., R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, P. Salvador, J.J. Dannenberg, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, A.G. Baboul, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, J.L. Andres, C. Gonzalez, M. Head-Gordon, E.S. Replogle, J.A. Pople, Gaussian 98, Revision A.11, Gaussian, Inc., Pittsburgh, PA, 2001. [4] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.J. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, J. Comput. Chem. 14 (1993) 1347–1363. [5] F. Pascale, C.M. Zicovich-Wilson, F. Lopez, B. Civalleri, R. Orlando, R. Dovesi, J. Comput. Chem. 25 (2004) 888–897. [6] E. Arroyabe, R. Kaindl, D.M. Többens, V. Kahlenberg, Inorg. Chem. 48 (2009) 11929–11934. [7] R. DeMichelis, B. Civalleri, M. Ferrabone, R. Dovesi, Int. J. Quant. Chem. 110 (2010) 406–415. [8] R. Kaindl, D.M. Többens, V. Kahlenberg, J. Raman Spectrosc. 42 (2011) 78–85.

[9] C. Adamo, V.J. Baronea, Chem. Phys. 110 (1999) 6158–6170. [10] E. Arroyabe, F. Prechtel, D.M. Többens, R. Kaindl, V. Kahlenberg, Eur. J. Mineral. (2011), doi:10 1127/0935-1221/2011/0022-2096. [11] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [12] P. Pulay, G. Fogarasi, G. Pongor, J.E. Boggs, A. Vargha, J. Am. Chem. Soc. 105 (24) (1983) 7037–7047. [13] A.P. Scott, L.J. Radom, Phys. Chem. 100 (1996) 16502–16513. [14] P. Borowski, A. Drzewiecka, M. Fernández-Gómez, M.P. Fernández-Liencres, T. ˜ Ruiz, Vib. Spectrosc. 52 (1) (2010) 16–21. Pena [15] R. Dovesi, V.R. Saunders, C. Roetti, R. Orlando, C.M. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N.M. Harrison, I.J. Bush, Ph. D’Arco, M. Llunell, CRYSTAL09 User’s Manual, University of Torino, Torino, 2009. [16] Y. Noel, M. Catti, Ph. D’Arco, R. Dovesi, Phys. Chem. Miner. 33 (2006) 383–393. [17] D.D. Johnson, Phys. Rev. B38 (1988) 12807–12813. [18] F. Pascale, C.M. Zicovich-Wilson, R. Orlando, C. Roetti, P. Ugliengo, R. Dovesi, J. Phys. Chem. B109 (2005) 6146–6152. [19] M.D. Towler, N.L. Allan, N.M. Harrison, V.R. Saunders, W.C. Mackrodt, E. Apra, Phys. Rev. B50 (1994) 5041–5054. [20] C.M. Zicovich-Wilson, F.J. Torres, F. Pascale, L. Valenzano, R. Orlando, R. Dovesi, J. Comput. Chem. 29 (13) (2008) 2268–2278. [21] R. Dovesi, C. Roetti, F. Fava, M. Prencipe, V.R. Saunders, Chem. Phys. 156 (1991) 11–19. [22] M. Causa, R. Dovesi, C. Pisani, C. Roetti, Phys. Rev. B33 (1986) 1308–1316. [23] CRYSTAL Basis Sets Library, 2006, http://www.crystal.unito.it/Basis Sets/Ptable.html [last accessed: 2010]. [24] L. Levien, C.T. Prewitt, D.J. Weidner, Am. Mineral. 65 (1980) 920–930. [25] J.P. Perdew, M. Ernzerhof, K. Burke, J. Chem. Phys. 105 (22) (1996) 9982–9985. [26] J.A. Pople, M. Head-Gordon, D.J. Fox, K. Raghavachari, L.A. Curtiss, J. Chem. Phys. 90 (10) (1989) 5622–5629. [27] D. Cremer, Z. He, J. Phys. Chem. 100 (1996) 6173–6188. [28] Y. He, D. Cremer, Mol. Phys. 98 (18) (2000) 1415–1432. ˜ [29] D.I. Bilc, R. Orlando, R. Shaltaf, G.-M. Rignanese, J. Íniguez, Ph. Ghosez, Phys. Rev. B77 (165107) (2008) 1–13. [30] A. Chopelas, Am. Mineral. 76 (1991) 1100–1109. [31] B. Kolesov, C. Geiger, Phys. Chem. Miner. 31 (2004) 142–154. [32] K. Iishi, Am. Mineral. 63 (1978) 1198–1208. [33] J.L. Servoin, B. Piriou, Phys. Status Solidi. B55 (1973) 677–686. [34] H. Suto, H. Sogawa, S. Tachibana, C. Koike, H. Karoji, A. Tsuchiyama, H. Chihara, K. Mizutani, J. Akedo, K. Ogiso, T. Fukui, S. Ohara, Mon. Not. Astron. R. Soc. 370 (2006) 1599–1606. [35] C.M. Zicovich-Wilson, F. Pascale, C. Roetti, V.R. Saunders, R. Orlando, R. Dovesi, J. Comput. Chem. 25 (2004) 1873–1881. [36] P.A.M. Dirac, Proc. Cambridge Philos. Soc. 26 (1930) 376–385. [37] S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200–1211. [38] J.P. Perdew, A. Ruzsinszky, G.T. Csonka, O.A. Vydrov, G.E. Scuseria, L.A. Constantin, X. Zhou, K. Burke, Phys. Rev. Lett. 100 (2008) 1–4, 136406. [39] R.K. Sato, P.F. McMillan, J. Phys. Chem. 91 (1987) 3494–3498. [40] F. Gervais, B. Piriou, Phys. Rev. B11 (1975) 3944–3950. [41] S. Ichikawa, J. Suda, T. Sato, Y.J. Suzuki, Raman Spectrosc. 34 (2003) 135–141. [42] N.G. Batalieva, Y.A. Pyatenko, J. Struct. Chem. 9 (5) (1968) 820. [43] S.P.S. Porto, R.S. Krishnan, J. Chem. Phys. 47 (3) (1967) 1009–1012.