ANNALS
OF PHYSICS:
67, 432-437 (1971)
Comments for Electronic
on the QUAD Band Structure
Scheme Calculations
G. GILAT* AND F. HERMAN IBM Research Laboratory, Monterey and Cottle Roads, San Jose, California Received December 28,197O
The QUAD and the Gilat-Raubenheimer schemes for calculating densities of states are discussed on a comparative basis. A semiquantitative analysis is drawn to enable such a comparison for tests of accuracy and convergence.
In a recent paper of this journal by Mueller, Garland, Cohen and Bennemann [l] (MGCB), a computational scheme called QUAD is introduced for calculating the electronic density of states g(E) in crystalline solids. In this procedure, one computes the n-th electronic band energy E,(k) at a relatively small number of mesh points k, in wave-vector space. This is accomplished by employing a suitable interpolation scheme involving neighboring k vectors. The main new feature of QUAD is that it includes not only the components of VE,(k,,), but also the six independent components of the symmetrized tensor I/a2E,/akiaki I(.These coefficients are utilized to generate in a random fashion (Monte Carlo) many more values of En(k) at points k in the vicinity of the various k, , thereby yielding a statistically improved density of states g(E). MGCB illustrate the QUAD procedure by computing g(E) for platinum (see Figs. 3 and 4 of Ref. [l]). In calculating this g(E), MGCB use 1876 matrix diagonalizations to generate about lo6 random points in k space. In their introduction, MGCB review earlier methods of calculating densities of states. In particular, reference is made to a method proposed by Gilat and Raubenheimer [2] (GR). In contrasting the QUAD scheme with the GR procedure, MGCB mention two main points of criticism of the GR method: (1) The GR approach uses only the linear terms of VE,(k) of the Taylor expansion. Therefore, the range of validity of the Taylor expansion about each point k is considerably smaller in comparison with the QUAD scheme which includes quadratic terms as well; (2) The GR method was originally applied to phonon spectra which usually are much less complex than electronic band spectra. Consequently, while the GR * Permanent Address: Department
of Physics, Technion, Haifa, Israel.
432
ELECTRONIC
BAND
STRUCTURE
433
CALCULATIONS
method is efficient in tackling lattice dynamical spectra, it may prove inadequate for more complex spectra which may include as many as 250 critical points. It is indeed true that in comparion to the QUAD scheme the GR method employs only the gradients VE,(k,J and neglects higher order terms, but there is a strong motivation for this. The GR method uses the linear expansion terms to obtain analytical expressions for the contributions to g(E) arising from a small volume centered at k, . This is, perhaps, the most powerful feature of the GR method, which is missing in the QUAD scheme. It may be pointed out that the GR method has recently been extended to the more general cases of the imaginary [3] and real [4, 51 parts of any spectral function in solids. It has also been extended to symmetries lower than cubic [6-81. Moreover, the GR method has already been used with considerable success in various calculations of the optical spectrum and electronic density of states of crystalline solids [9-121. We now proceed to make a more critical comparison of the two methods. MGCB argue that by including quadratic terms it is possible to obtain better approximations to E(k) over a wider range at k, and thus increase substantially the domain of the expansion. This would lead, they claim, to computer time saving by reducing the number of necessary diagonalizations of the Hamiltonian matrix. MGCB are probably right, in principle, in claiming that the second-order approximation is better than the first-order one. The only trouble is that it is also by far less convenient. In going beyond the linear approximation, MGCB sacrifice the analytical nature of the calculation and resort to the much inferior Monte Carlo technique. The question is whether this sacrifice is really justified. The best answer to this would be to perform the same calculation with both techniques and compare results. Unfortunately, such a comparative computation has not been carried out yet, so that it is not easy to draw a thorough and critical comparison. Nevertheless, MGCB make a few off-hand estimates according to which the range of reliability of each k, point in the QUAD scheme is 5 times larger than that of the GR scheme. They conclude, therefore, that it would take 375 000 matrix diagonalizations for the GR method to achieve the same accuracy for g(E)as the QUAD does with only 1876 diagonalizations. We consider these estimates as highly questionable. In order to obtain more realistic estimates for comparison purposes, we propose the following semiquantitative analysis. Letf(x) be a function given for a domain 0 < x < 1 and for which we want to make an estimate over a range a < 1, i.e., 0 < 1x - x0 1 < a. We perform the Taylor expansion of f(x) and, for comparison, we consider the zero-th, first- and second-order terms, respectively:
f(x) = .h + ... f(x) =h + hf’ + ...
(0 W, (1 St),
f(x)
(2
=fo
+ hf’ + h2f”/2 + ..+
n4,
434
GILAT
AND
HERMAN
where f0 = f(x,J and h = x - x,, . Next, we specify an “error” 6 within which f(x) is to be approximated, respectively, for the various expansions. We let a,, , a, and a2 denote the respective ranges within which this is to be accomplished. We also introduce the linear mesh numbers n, = l/ai . In analyzing this example we notice that the leading term in 6 for the 0-th order expansion is the first-order term, namely hf’. We can, therefore, estimate [13] a, , or rather n, by 1f ‘/S 1. By similar reasoning, we obtain the following estimates for ni : f (1l/2 f” l/3 no = and nl=s I I , (2) n2=G I I * Next we look for a possible link between the various ni . This can be found by noting, that in the problems usually treated,f(x) is periodic in x, and can, therefore, be Fourier expanded over the domain 0 < x < 1. We assume therefore, that f(x)= C+olcos7rx+psin7~x++** where for simplicity only first-order components are considered. It is clear that without knowing 01and /I it is impossible to draw quantitative conclusions. Nevertheless, because of the periodic nature of f(x), one can estimate the mean value of ( f(n-l)/f(n)) over the entire range of 0 < x < 1, and it is l/n. Substituting for this in Eq. (2), we obtain n2 : n, : n, = ($ n,)1’3 : (5 n,)l”
: (no)‘.
(3)
On a three-dimensional grid, the number of mesh points mi which are equivalent to the ni , respectively, are related by m2 :m 1 : m. s n, 1 : #Z : no’,
(4)
where all the constant coefficients are left out for simplicity. Since this result is based on very general and simplified grounds, we must caution the reader not to accept it too literally. Nevertheless, regardless of the dubious numerical factors in Eq. (3), the powers of l/3, l/2, and 1 are significant and may serve as a basis for an objective comparison of the methods. It is quite obvious that m, is associated with the number of mesh points required for the root sampling (histogram) method, whereas ml and m2 are associated with the GR and QUAD schemes, respectively. Using Eq. (4) we now demonstrate to what extent the Monte Carlo technique is inferior to the analytical method. In the Pt computation [l], MGCB performed about 2000 matrix diagonalizations to obtain a certain degree of accuracy in g(E). Using Eq. (4), where m2 N 2000, one obtains for m,, N 10s-lO1o mesh points. Nevertheless, MGCB were able to generate only 10Gmesh points and the reason for this is their use of the Monte Carlo technique. In other words, relative to analytical integration, the efficiency of the Monte Carlo
ELECTRONIC
BAND
STRUCTURE
CALCULATIONS
435
technique is about 1 : 1000 or less. In comparison, in order to attain an accuracy equivalent to log root-sampling points, the GR method would require about m, = 6 x lo4 points. However, for an accuracy of lo6 root-sampling points, the GR method requires only about 2 x lo3 matrix diagonalizations, very similar to the QUAD scheme. The conclusion of this discussion is that anything that might be gained by including quadratic terms is liable to be lost again by the subsequent use of Monte Carlo techniques. It is interesting to point out that, in a recent paper, Dalton [14] tests for convergence of the GR method for the cases Fe, Cu and Ni. Although he does not use more than about m, = 400 mesh points, his results suggest that the quality of his g(E) graphs may resemble that of MGCB’s for an equal number of mesh points. The second point of criticism of MGCB of the GR method concerns the accuracy near critical points, where V&(k) is quite small. According to MGCB, this may seriously affect computations of the density of states for cases where one may expect as many as 250 Van-Hove singularities. We do not intend to discuss this point in detail here, but it is obvious that if there is a problem here, QUAD is not likely to be its solution. This can easily be inferred by looking at Figs. 3 and 4 of Ref. [l]. It is very hard to envisage anything like 250 critical points in these plots of g(E). Their number is more in the vicinity of 20 to 25, and it is not easy to give an exact figure because many of the critical points are blurred by the somewhat crude resolution of these graphs. Thus, their criticism here is unconvincing. In evaluating the QUAD scheme, it seems to us that it is severely handicapped by using statistical methods rather than analytical integration. Unfortunately, the reason for this is the inclusion of quadratic expansion terms which otherwise should improve matters. A good way to show the seriousness of this handicap is to consider the computation of properties such as Ed, or more generally, the principal (real) value of linear response (Green) functions. In their paper, MGCB admit that this particular task is beyond the reach of QUAD. The source of their difficulty is the singular character of the integrand which, if performed by statistical methods such as the Monte Carlo technique, is bound to bring in a large amount of statistical error. This difficulty is readily overcome by using analytical integration, as was recently shown by Gilat and Bohlin [4, 151. Another possible weak point in the QUAD scheme is their method of interpolation. The danger of this arises from the possibility of wrongly identifying the various bands for the purposes of performing the quadratic expansion. In the QUAD scheme, the band identification is made according to En-, < En < En+, . In regions, where two bands rapidly approach each other, this may lead to errors in estimating the gradients. It is hard to estimate the significance of this error, but it is believed to be relatively less important. In conclusion, the QUAD scheme, although it has many remarkable features, does not seem to be the remedy for the shortcomings that its authors ascribe to
436
GILAT AND HERMAN
earlier
be
k,
main
about E(k,), but the best features
to generate
each
Its
techniques.
expansion combining
many
more
feature
is the
the price of the values
inclusion
of E,(k)
and
in the Taylor An elegant way of methods would, perhaps,
of quadratic
paid for this is too GR and the QUAD V&(k)
terms
high.
on a reguZur
mesh
of k about
, where
E,(k) = K&,) + VWo) . (k - ko) + (k - k&2 - Q - (k - kc,), VE,@) = V&&J + Q * (k - ko),
(5) (6)
and Q is the symmetrized tensor I/ a2En/aki ak, jj . These additional mesh points k can then serve as centers of rectangular cells which properly dovetail into each other to fill the entire volume of integration [16]. The GR method can then be applied
individually
to each
of these
cells
to perform
the
integrals
analytically.
Notes added in proof. (1) In a computation performed after this paper was submitted for publication, Cooke and Wood (to be published) perform detailed comparisons of the GR, the QUAD and the “combined” schemes [16]. The comparisons of Cooke and Wood strongly support the arguments of the present article. (2) It was brought to our attention that Brust [17] had employed a method very similar to the QUAD scheme in his calculation of photoemission from Si.
ACKNOWLEDGMENT It is a pleasure (G. G.) wishes enabling him to Trade Research kind hospitality
to thank Dr. Norris W. Dalton for many illuminating discussions. One of us to express his gratitude to IBM World Trade Corporation and IBM-Israel for spend a period of time at the IBM San Jose Research Laboratory as a World Fellow. He is also grateful to the IBM San Jose Research Laboratory for the he enjoyed.
REFERENCES 1. F. M. MUELLER, J. W. GARLAND, M. H. COHEN, AND K. H. BENNEMANN, Ann. Phys. See also F. M. MUELLER in “Computational Methods of Band Theory” p. 305, International Conference, 1970 (P. M. Marcus, J. F. Janak and A. R. Williams eds.) Plenum Press, New York, 1971. 2. G. GILAT AND L. J. RA~ENHEIMER, Phys. Rev. 144 (1966), 390. 3. G. GILAT AND Z. KAM, Phys. Rev. L&t. 22 (1969), 715. 4. G. GILAT AND L. BOHLIN, Solid State Comm. 7 (1969), 1727. 5. N. W. DALTON, Solid State Comm., to appear. 6. L. J. RAUBE NHEIM@R AND G. GILAT, Phys. Rev. 157 (1967), 586. 7. Z. KAM AND G. GILAT, Phys. Rev. 175 (1968), 1156. 8. E. FINKMAN, Z. KAM, E. COHEN AND G. GILAT, to be published. 9. C. W. HIGGINBOTHAM, F. H. POLLAK AND M. CARDONA, Solid State Comm. 5 (1967), 916.
ELECTRONIC
BAND
STRUCTURE
CALCULATIONS
437
10. L. R. SARAVIA ANO D. BRUST, Phys. Rev. 170 (1968), 11. F. HERMAN, R. K. KORTIJM, C. D. KUGLIN, AND J.
683. L. SHAY, in “II-VI Semi-Conducting Conference, 1967, (D. G. Thomas, Ed.) Benjamin, New
Compounds,” p. 503, International York, 1967. 12. J. F. JANAK, D. E. EASTMAN AND A. R. WILLIAMS, Solid State Comm. 8 (1970), 271. See also J. F. JANAK, in “Computational Methods of Band Theory” p. 323, International Conference, 1970 (P. M. Marcus, J. F. Janak and A. R. Williams eds.) Plenum Press, New York, 1971. 13. The various ranges ai are perhaps somewhat overestimated, but since this is meant to be only a guide of thought, this error should not be taken seriously. 14. N. W. DALTON, J. Phys. C. 3 (1970), 1912. 1.5. N. W. DALTON AND G. GILAT, Bull. Amer. Phys. Sot. 15 (1970), 1635. 16. A similar approach is currently being used by J. F. COOKE AND R. F. WOOD (to be published) see also [12]. 17. D. BRUST, Phys. Rev. 139 (1965), A489.