CRYSTAL GROWTH & DESIGN
Quantum-Mechanical and Thermodynamical Study on the (100) Face of LiF Crystals
2009 VOL. 9, NO. 1 404–408
Marco Rubbo,* Marco Bruno, and Mauro Prencipe Dipartimento di Scienze Mineralogiche e Petrologiche, UniVersita` degli Studi di Torino, Via Valperga Caluso 35, I-10125 Torino, Italy ReceiVed June 12, 2008; ReVised Manuscript ReceiVed September 23, 2008
ABSTRACT: In this work we determine the Helmholtz free energy (F) of the (100) face of LiF. It accounts for the athermal work of separation at 0 K, the vibrational and zero point surface energy, and the surface entropy. The thermal contributions to F are calculated via the vibrational partition functions obtained from the frequency spectrum of the infinite crystal and of slabs limited by (100) faces. The vibrational frequencies and modes are calculated using the CRYSTAL06 software for quantum-mechanical ab initio calculations. The relevant surface properties are the thermodynamic excess functions of a macroscopic slab, limited by the face of interest, in comparison with the bulk crystal. Strengths and weaknesses of the method used are discussed. Introduction We are here concerned with the determination of both the stable faces making the equilibrium habit of a crystal and of their surface structure. Rather than the surface free energy at the temperature of interest, in order to determine the stability of a crystal face, the athermal work of separation at T ) 0 K is usually evaluated. In this way, one can overlook the morphological importance of faces experiencing both relaxation and reconstruction: as a case study we refer to halite (NaCl) and, in particular, to a recent work on this mineral1 and to references therein. Determining the stability of a face leads one to consider several different surface structures which are potentially stable, thus allowing for (i) the rationalization of the formation of epitactic layers on crystallographic and energetic grounds, (ii) the identification of surface site candidates for foreign adsorption, (iii) the prediction of a possible change of the character of a face from kinked (or stepped) to stable,2,3 and (iv) the analysis of interfacial phenomena on a sound starting point, in general. For these reasons, the less stable crystal faces are also clearly worth of a characterization. The thermal contribution to the energy of a crystal face can now be calculated using powerful computer codes: in this work we aim to test the potential of CRYSTAL064-6 and its weight in terms of computer time spent to evaluate the thermodynamic functions of a crystal face by a quantum-mechanical ab initio calculation. Being interested in the properties of alkali halides, in particular in their equilibrium and growth habit, we started by assessing the contribution of the static and vibrational energies, and entropy to the free energy of the (100) face of LiF (space group Fm3jm, Z ) 4). As a preliminary work, in order to obtain a good compromise between accuracy and computing time, we performed tests on the bulk crystal by using different Hamiltonians and basis set functions. With the choice of the computational parameters we made (see Appendix), we obtained for the bulk crystal a value of the optimized conventional cell parameter equal to 4.01 Å, which is in excellent agreement with the experimental one (a * Corresponding author. E-mail:
[email protected], tel. +390116705127, fax: +390116705128.
Figure 1. Variation of the area of the 2D (100) cell (in Å2) with the number of layers in the slab. Diamonds correspond to calculated cell area. The continuous line corresponds to the least-squares fit: A(n) ) 8.0309 - 1.9331(n-1) + 2.6372(n-2) - 3.0921 (n-3) [Å2] (see main text for details).
) 3.99 Å7) at 298.15 K. See also the work of Prencipe et al.8 for discussions of correlation effects in alkali halides. Results and Discussion The CRYSTAL06 code allows for the calculation of surfaces thermodynamic functions by comparing the bulk properties with those of a two-dimensional periodic structure (henceforth termed slab) of given thickness, parallel to the crystal plane of interest. This study deals with slabs limited by (100) planes. The slab unit cell contains a number of LiF formula units equal to the number of layers in the slab. The thinnest slab considered comprises four layers. For the lack of translational symmetry in the direction [001], after optimization this slab has a thickness of 6.0806 Å and 2D mesh parameters a ) b ) 2.7684 Å, R ) 90, the vectors a and b being directed along the diagonals of the 2D mesh of the conventional cell. The neighbor Li-F distances vary slightly from layer to layer, after optimization. The thinner the slab the higher the shrinking of the 2D mesh and the structural relaxation; as the thickness of the slab increases, the 2D mesh tends to reach the bulk size. The surface area of the 2D mesh as a function of the number of layers are represented in Figure 1.
10.1021/cg800615y CCC: $40.75 2009 American Chemical Society Published on Web 12/04/2008
Study on the (100) Face of LiF Crystals
Crystal Growth & Design, Vol. 9, No. 1, 2009 405
Figure 3. Average distance between subsequent layers, from the top one (i ) 1) to the middle of the slab 28 layers thick.
inward and F anions are shifted outward, but in the underneath layers ∆i(100) converges rapidly to zero (Figure 2b). Within the slab of 28 layers, the average distance d i(100) between layers, after some oscillations of rapidly decreasing amplitude, reaches an average value of 2.0104 Å for the layers close to the center of the slab; the latter value is higher than the optimized value in the LiF bulk (2.0035 Å) by about 3.4% as shown in Figure 3. The variations of ∆i(100) and d i(100) are qualitatively similar to those calculated in the case of NaCl by Bruno et al.1 Surface Energy. The specific surface energy, Es (erg/cm2), was calculated with the following relation:10 Figure 2. (a) The relaxed structure of the (100) slab viewed along the direction [001]. (b) The variation of the absolute value of rumpling from the slab surface to the bulk. The slab comprises 28 layers.
Geometry optimizations (lattice parameters and atomic coordinates) of the (100) slabs up to a thickness of 28 layers were performed; at the maximum thickness considered, the bulk-like properties in the central zone of the slab are closely approached. Listing of lattice parameters and atomic coordinates of the slabs are available on request. The shrinking of the 2D mesh is accompanied by a variation of the interlayer distance along the [001] direction: the Li and F atoms that make the outer layers define parallel planes at different heights with respect to the slab center; as a result, the mean equidistance is changing along and with the thickness of the slab. However, for sufficiently thick slabs, the outer layer configuration no longer changes as the number of layers is increased. To describe the thickest (100) slab structure we use the following parameters:9
1 ∆i(100) ) (xiLi - xiF) 2 1 Li F + xi+1 - xiLi - xiF) di(100) ) (xi+1 2
Es ) lim Ee(n) ) lim nf∞
m
of the kind: P ) ∑ cin1-i i , where ci are constants. i)1 In Figure 4 the surface energy is plotted as a function of the number of layers. The dots are the calculated energies of the relaxed surfaces and the line is the continuous representation:
Ee(n) ) 350.9 -
where xiLi and xiF are the absolute coordinates along [001] of Li and F, respectively, in the ith layer (for i ) 1, . . ., n); d i(100) is the average distance to the next underlying layer, and ∆i(100) is the average rumpling of the ith layer, the latter corresponding to the vertical displacements, of equal magnitude but opposite sign, of the cations and anions, off the geometrical average of the layer. A positive rumpling is connected to an outward movement of the cations, and a negative one to an inward motion. In Figure 2a the optimized (100) surface structure is schematically reported. At the surface the Li cations are shifted
(3)
where E(n)slab is the energy of an n-layer cell comprising n formula units; A(n) is the area of the primitive unit cell, as defined before; Ebulk is the energy of a formula unit in a bulky crystal; Ee(n) the energy per unit area required to form an n-layers slab. As more layers are added in the calculation (nf∞), Ee(n) will converge to the surface energy per unit area, Es. We estimated the thermodynamic properties by fitting the calculated quantities with functions of the number of layers (n) in the slabs: this assumes (a) the existence of continuous functions taking on the correct values of the property of the assembly of layers (an integer number) making the slabs, and (b) that these functions can be represented by truncated series
(1) (2)
nf∞
E(n)slab - nEbulk 2A(n)
388.7 1549.7 2574.7 + n n2 n3
(4)
The surface energy of a thick slab so calculated coincides with the athermal surface energy at 0 K and zero pressure. The surface energy of the 28 layers slab is 338.9 erg/cm2, whereas the surface energy extrapolated for a slab of macroscopic thickness is Es ) 350.9 erg/cm2. These values can be overestimated as a consequence of the basis set superposition error.6 We performed some tests using the technique of the “ghost atoms”:6 that is a layer of atoms is added, deleting the nuclear and electronic charges, but leaving the basis set centered at the atomic positions. We estimate that the error on the calculated energy is approximately +5%. Such error would not affect the
406 Crystal Growth & Design, Vol. 9, No. 1, 2009
Figure 4. Surface energy [erg/cm2] of the (100) slab as a function of thickness.
calculation of the vibrational spectrum, as the latter is based on the derivatives of the energy. By means of quantitative experiments, Gilman11 measured the reversible work to cleave LiF-forming (100) surfaces, at 77 K. He obtained a value of 340 erg/cm2, which represents the free energy of the crystal related to the formation of a unit area of free surface. It includes therefore the energy of vibration of the newly formed surface atoms, together with the energy due to the relaxation of their positions with respect to the bulk.12 In order to be able to compare the experimental and measured values we have to include the thermal energy and entropy whose contributions are calculated in the following. Frequency Calculation. CRYSTAL0613,14 computes numerically the second derivatives of the energy using analytical first derivatives. The calculation is done by displacing the irreducible atoms only. The full mass-weighted Hessian matrix, in Cartesian coordinates, is then generated and diagonalized to obtain eigenvalues and normal modes. The frequencies must be computed when the crystalline structure is at a minimum on the potential energy surface. Therefore, harmonic frequency calculations have been performed on the optimized (100) slabs, which were built on the basis of the conventional cell vectors; five slabs of increasing thickness, given as multiple of the conventional cell parameter in the [001] direction, were considered. The 2D slabs are made of 4, 8, 12, 16, and 20 layers. To compute the excess thermodynamic functions, the frequency calculations were also performed on bulk crystals built by supercells with the same number of layers in the [001] direction as the slabs. Such crystals have a P4/mmm space group; the space group of the 2D slabs is P4mm. In general, the decreased number of symmetry operators (and the consequent partial degeneracy removal), together with the (100) truncation effect (surface), modify the frequency spectra at each point of the Brillouin zones of the 3D crystal and 2D slabs. The comparison of the frequency distributions is meaningful provided that slabs and bulky crystal have the same supercell. In Figure 5 the slab is 16 layers thick and the size of the bulk supercell is 4a × a × a where a is the parameter of the conventional fcc cell. In comparison with the bulk, the (100) slabs exhibit a distribution richer in modes at low frequency, and a higher dispersion. From the analysis of the eigenvectors of the massweighted Hessian matrix (normal modes, which are available upon request), it appears that long wavelength modes (short |k|) have similar patterns in the slab and in the bulk: the atoms move without alteration of the distances within each reticular plane
Rubbo et al.
perpendicular to the wave vector, while the planes move in phase. As the wavelength is decreased, modes get excited where atoms of the same charge in neighbor rows move out of phase. The highest frequency transverse modes can be described as stretching of the neighbor Li-Li distances accompanied by a minor deformation of the F-F distances, producing an overall distortion of the basic octopole. A distortion of the basic octopole also occurs in the longitudinal highest frequency modes having a wave vector parallel to [001]. Slabs show typical surface modes characterized by the abrupt decrease of the amplitude of the oscillations from the surface toward the bulk of the slab. This is less evident in the thinner slabs, but it is clearer when slabs have 12 or more layers. More precisely, four surface modes of E symmetry appear in the interval 304 and 322 cm-1 when slabs consist of 20, 16, and 12 layers, respectively. For such E modes, the oscillations of the atoms at the surface are parallel to the plane (100) and occur in the directions [100] and [010]; the Li atoms oscillate with higher amplitude and in antiphase in respect to the F atoms. Two surface modes, with A and B symmetry, occur at 188-189 (12 layers), 191-192 (16 layers), 193-194 cm-1 (20 layers); in these cases the F ions oscillate in a direction perpendicular to (100), out of phase with both the neighbor F ions, in the same plane, and with neighbor Li ions in the direction [100]. Evaluation of Thermodynamic Properties. The frequency spectrum of a bulk crystal and that of a slab limited by the faces of interest determine different vibrational partition functions and thermodynamic properties. The difference between the values for the slabs and for the bulk crystal, of a given thermodynamic property, become constant when the slab reaches a macroscopic size: it corresponds to the desired value of the surface property. Therefore, to evaluate this limiting value, the frequencies and vibrational modes have been calculated, at the Γ (k ) 0) point, for slabs with increasing thickness in the [001] direction. In parallel the thermodynamic properties of bulk crystals, built by supercells as thick in the direction [001] as the slabs, were calculated. The lattice parameters of the supercells are: a ) b ) 4.01 Å, c ) L × 4.01 Å, R ) β ) γ ) 90°; the reciprocal space group is P4/mmm* (the asterisk indicates an entity in reciprocal space). The vibrational modes and frequencies of the bulk crystal generated by such supercells correspond either to modes at the Γ point or at k * 0 in the direction ∆, i.e. [001]*. For the thickness chosen (corresponding to L ) 1, 2, 3, 4, 5), we obtained frequency values at points k ) [0, 0, (2πm)/(La)], m ) 0, 1, . . ., L - 1; the resultant dispersion curves are shown in Figure 6, where 24 frequency values for each k point are plotted. At the Γ point, the IR active frequencies corresponding to the longitudinal modes experience the so-called TO/LO shift (TO and LO stand for transverse optical and longitudinal optical mode, respectively). By knowing the optical dielectric constant15 of LiF (1.96), such a shift has been calculated with CRYSTAL06 4,6 and amounts to 361 cm-1 for the unique IR active mode of F symmetry. To correctly compare the properties of the slabs with those of the bulk crystals, for the latter one we have shifted by such an amount only the frequency of the IR active LO mode, at the Γ point, which has the wave vector parallel to [001]*. Indeed, in the slabs, the lack of translational symmetry along [001], and the consequent removal of the Born-Von Karman boundary conditions, results in frequencies already automatically corrected for the LO/TO shift for those modes having wave vectors parallel to [001]*; by contrast, the frequencies corresponding to all of the other modes should be corrected, as in the case of the bulk crystal, provided the dielectric constants of
Study on the (100) Face of LiF Crystals
Crystal Growth & Design, Vol. 9, No. 1, 2009 407
Figure 5. Number of modes vs frequency in a slab (left) and in the bulk (see text for detail). The resolution of the histograms is 7.5 cm-1.
Figure 6. Dispersion curves as a function of the reduced wave vector in the ∆ direction. The highest frequency at Γ, corresponds to an LO mode.
Figure 7. Stars are calculated excess entropy of slabs of increasing number of layers, at the temperature T ) 77 K. The continuous line is a least-squares fit given in the text.
the slabs are known. However, since the latter quantities are unknown, for the sake of slab-bulk comparison, the frequencies of the latter modes were uncorrected both in the slab and in the bulk, thus assuming an error cancelation when excess properties, depending upon the difference of the frequencies in the slabs and bulk crystals, are calculated. Entropy was calculated for the bulk crystals and slabs from the partition functions associated to the respective normal modes:
entropy per unit area, Ss. It is worth observing that the increase of the surface area with the slab thickness contributes to the decrease of both surface energy and entropy. We calculated the surface entropy at 77 K for the sake of comparison with the experimental value of the surface free energy by Gilman.11 For the thicker, 20-layer slab it amounts to 0.1081 erg/(cm2K). The continuous line in Figure 7 corresponds to Se(n) ) 0.1845 - (1.9467)/(n) +(9.034)/(n2) - (15.8571)/(n3) erg/(cm2 K); the extrapolated value of the surface entropy of a (100) crystal face is SS ) 0.1845 erg/(cm2 K) at 77 K. We found that the variation with temperature of the harmonic excess entropy of the 20-layer slab depends weakly on temperature. It increases up to ≈ 0.12 erg/cm2 K at 120 K and then decreases toward ≈0.10 erg/cm2 K at higher temperatures; however, with increasing temperature anharmonic effects should be considered. Vibrational Energy. To compute the surface free energy for the slabs, Fe(n) (erg/cm2), and to obtain the extrapolated value for a macroscopic face, Fs, we first calculated the zero point and thermal vibrational energy. The thermal contribution to Fe(n) is defined in eq 8:
S ) -kB
[
(
ν) ∑ hω(k, ν)ln 1 - exp - hω(k, kBT k,ν
)]
+
∑
1 hω(k, ν)n(ω, T) T k,ν
(5)
where the number of phonons in the ν branch with wave vector k follows the Bose-Einstein distribution:
[ (
n(ω, T) ) exp
) ]
hω(k, ν) -1 kBT
-1
(6)
In the two last relations h and kB are the Planck and Boltzmann constant, respectively, and ω is the angular frequency of the wave (k,ν). Finally, the (100) specific surface entropy, Ss (erg/cm2 K), was calculated by subtracting from the entropy of slabs of different thicknesses, S(n)slab, that of an equal number of atoms in the bulk, S(n)bulk (eq 7).
S(n)slab - S(n)bulk Ss ) lim Se(n) ) lim nf∞ nf∞ 2A(n)
Fth_e(n) )
F(n)th_slab - F(n)th_bulk 2A(n)
(8)
F(n) associated to slab or bulk was calculated by the relation:
(7)
In analogy with the surface energy calculation, as more layers are added in eq 7 (nf∞), Se(n) will converge to the surface
F(n)th_slab/bulk ) E(n)th_slab/bulk - TS(n)slab/bulk
(9)
where E(n)th_slab/bulk includes the zero point energy according to eq 10:
408 Crystal Growth & Design, Vol. 9, No. 1, 2009
E0_slab/bulk + Eth_slab/bulk )
Rubbo et al.
∑ hω(k, ν){ 21 + n[ω(k, ν), T]} k,ν
(10) Then the excess E0_e + Eth_e is the difference between the vibrational energies (obtained by summing over the respective vectors k and branches ν, eq 10) of n-layer slabs and bulk crystals having the appropriate supercell. We obtained for the 20-layer slab, at 77 K, the following values: zero point energy E0_e ) 2.4029 erg/cm2, thermal energy contribution Eth_e ) 4.0706 erg/cm2, and free energy of vibration Fth_e ) E0_e + Eth_e - TS_e ) -1.8492 erg/cm2. At 77 K,
Fth_e(n) ) -5.6510 +
108.999 646.4522 1837.7134 + n n2 n3 (11)
so the extrapolated value of the surface vibrational free energy is Fth_s ) -5.6510 erg/cm2, for the (100) face. Finally, adding the vibrational free energy Fth_s to the athermal surface energy previously calculated at 0 K, the Helmholtz free surface energy of (100) LiF amounts to Fs ) 345.2 erg/cm2 at 77 K in excellent agreement with the measured11 340 erg/cm2. With increasing temperature the weight of the entropy term becomes more important, determining a value of Fth_s equal to -47.1225 erg/cm2 for a macroscopic slab at 298.15 K. Conclusions From ab initio quantum-mechanical calculations, we have obtained a value of the surface free energy of the (100) face of LiF in a quite satisfactory agreement with the experimental datum obtained by Gilman.11 This result is an indication of the strength of the algorithms used6 and implemented in CRYSTAL06. From a fundamental point of view, neglecting thermal contribution introduces an error of ≈ 2% on the Helmholtz surface energy at 77 K, but the error increases to ≈ 13% at 298 K. Moreover, one is limited to a static picture of the crystal face and we think that for polar faces such as the (111) of LiF, the analysis of the normal modes of vibration is even more important for it allows assessment of the stability of the reconstructed surface. At variance with other methods (see for instance Parlinski16 and references therein), the approach to the study of surfaces of CRYSTAL064 is peculiar in that the relevant properties of 3D and 2D systems are treated on an equal footing which however demands a great amount of computing time. To cope with this constraint, we have successfully extrapolated the relevant thermodynamic properties of the face from the calculated values for the slabs. Appendix Computational Details. Bulk and slab geometry optimizations were performed by means of the ab initio Hartree-Fock LCAO SCF CRYSTAL06 computer program,4 for the study of periodic systems.5 The multielectronic wave function is constructed as an antisymmetrized product (Slater determinant) of monoelectronic crystalline orbitals (CO) which are linear combinations of local functions (to be indicated as AO’s) centered on each atom of the crystal. In turn, AO’s (basis set) are linear combinations of Gaussian-type functions (GTF, the product of a Gaussian times a real solid spherical harmonic to give s-, p-, and d-type AO’s).
In the present case, Li is described with a 6-11G* basis.17 It consists of six contracted GTF’s for the description of the core s shell, and two uncontracted sp and one d functions for the valence shell. For F, the basis sets used is 7-311G*18 consisting of seven contracted GTF’s for the description of the core s shell, and a contraction of three and two uncontracted sp’s, and one d function for the valence shell. All of the exponents of the valence sp and d GTF’s have been variationally optimized. The calculations were performed at the Hartree-Fock level. The five parameters (ITOL1, ITOL2, ITOL3, ITOL4, and ITOL5) controlling the accuracy of the calculation of Coulomb and exchange integrals14 were set to 10-8 (ITOL1 to ITOL4) and 10-16 (ITOL5). The diagonalization of the Fock matrix was performed at 10 k points in the reciprocal space (Monkhrost net19) by setting the shrinking factor IS14 to 8. For both geometry and frequency calculations, the SCF convergence threshold on total energy (TOLDEE in CRYSTAL06) has been set to 10-10 bohr. Geometry optimization is achieved when each of the components of the gradient and of the displacement (TOLDEG, TOLDEX parameters in CRYSTAL06) are smaller than 4 × 10-5 hartree bohr-1 and 10-5 bohr, respectively. Acknowledgment. The authors thank Professor Dino Aquilano for discussions and reading the manuscript.
References (1) Bruno, M.; Aquilano, D.; Pastero, L.; Prencipe, M. The structures and surface energies of (100) and octopolar (111) faces of halite (NaCl): an ab initio quantum-mechanical and thermodynamical study. Cryst. Growth Des. 2008, 8, 2163–2170. (2) Pastero, L.; Costa, E.; Bruno, M.; Rubbo, M.; Sgualdino, G.; Aquilano, D. Cryst. Growth Des. 2004, 4, 485–490. (3) Pastero, L.; Aquilano, D.; Costa, E.; Rubbo, M. J. Cryst. Growth 2005, 275, e1625–e1630. (4) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; Zicovich-Wilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M. CRYSTAL06 User’s Manual; University of Torino: Torino, 2006. (5) Pisani, C.; Dovesi, R.; Roetti, C. Hartree-Fock ab-initio treatment of crystalline systems, Lecture Notes in Chemistry; Springer: Berlin, 1988. (6) Basic theory, algorithm, and applications can be found at http:// www.crystal.unito.it. (7) The Handbook of Chemistry and Physics, 67th ed.; West, R. C., Ed.; CRC Press: Boca Raton, FL, 1986. (8) Prencipe, M.; Zupan, A.; Dovesi, R.; Apra`, E.; Saunders, V. R. Phys. ReV. 1995, B51, 3391. (9) Vogt, J.; Weiss, H. Surf. Sci. 2001, 491, 155–168. (10) Dovesi, R.; Civalleri, B.; Orlando, R.; Roetti, C.; Saunders, V. R. Ab initio quantum simulation in solid state chemistry. In ReViews in Computational Chemistry; Lipkowitz, B. K.; Larter, R.; Cundari, T. R., Eds.; John Wiley & Sons, Inc.: New York, 2005; Vol. 21, pp 1-125. (11) Gilman, J. J. J. Appl. Phys. 1960, 31, 2208–2218. (12) Mutaftschiev, B. The Atomistic Nature of Crystal Growth. In Springer Series in Material Science; Springer, Berlin, 2001; p 61. (13) Pascale, F.; Zicovich-Wilson, C. M.; Lopez Gejo, F.; Civalleri, B.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 888–897. (14) Zicovich-Wilson, C. M.; Pascale, F.; Roetti, C.; Saunders, V. R.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 1873–1881. (15) Physics of II-VI and I-VII Compounds; Hellwege, K.-H.; Madelung, ¨ rnstein, New Series, Group III, Springer: New O., Eds.;Landolt-BO York, 1987; Part b, Vol. 17. (16) Parlinski, K. Phys. ReV. 2006, B74, 184309. (17) Ojama, L.; Hermansson, K.; Pisani, C.; Causa`, M.; Roetti, C. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1994, 50, 268– 279. (18) Nada, R.; Catlow, C. R. A.; Pisani, C.; Orlando, R. Simul. Mater. Sci. Eng. 1993, 1, 165–187. (19) Monkhrost, H. J.; Pack, J. D. Phys. ReV. B 1976, 13, 5188–5192.
CG800615Y
