VOLUME 69, NUMBER 19 PH YSICAL REVIEW LETTERS 9 NOVEMBER 1992 Ab Initio Calculation of Force Constants and Full Phonon Dispersions Siqing Wei and M. Y. Chou School of Physics, Georgia Institute of Technology, Atlanta, Georgia 30332 043-0 (Received 2 June 1992) We present a method to calculate the full phonon spectrum using the local-density approximation and Hellmann-Feynman forces. By a limited number of supercell calculations of the planar force constants, the interatomic force constant matrices are determined. One can then construct the dynamical matrix for any arbitrary wave vector in the Brillouin zone. We describe in detail the procedure for elements in the diamond structure and derive the phonon dispersion curves for Si. The anharmonic eff'ects can also be studied by the present method. PACS numbers: 63.20. Dj The energy dispersion of phonons has long been a focus of interest because it provides rich information about the dynamical properties of the material [1]. In particular, it is an essential input in the calculation of heat capacities, thermal expansion coefficients, electron-phonon interac- tions, etc. , for the crystal. Attempts have also been made to study the vibrational properties of complex systems, such as semiconductor alloys and superlattices, based on information about the force constants of the constituent materials [2, 3]. Ab initio calculations of the force con- stant matrices and phonon spectra will not only give an accurate database for these applications but also provide stringent tests of various empirical models [4,5]. There are basically two methods in use for the calcula- tions of phonon frequencies within the framework of the local-density approximation (LDA): (1) the linear re- sponse theory with dielectric screening in which atomic displacements are treated as perturbations [6] and (2) the "direct" approach which calculates the total energy of the distorted system or the Hellmann-Feynman forces [7] on the atoms using the supercell method. Each of these methods has its advantages and drawbacks. In the first approach, the response to perturbations was calculated in the past by inverting the dielectric matrix [6] which is computationally cumbersome and restrictive. Recently, new schemes have been proposed to obtain the linear response either by iterating up to self-consistency [8] or by solving an integral equation [9] for the change in the electron density. These methods can handle perturba- tions of arbitrary wave vectors, yet only linear effects are considered. On the other hand, the direct approach which considers periodic distortions using supercells is computationally straightforward. It handles the perturbed and unper- turbed systems on the same footing under the frozen- phonon approximation, and allows one to study, in princi- ple, both linear and nonlinear effects. Phonon frequencies for isolated symmetry points can be easily calculated [10] by the pseudopotential LDA method. However, the su- percell size increases rapidly as the symmetry decreases. Only the dispersions along a few high-symmetry direc- tions have been reported [11-14] using this direct ap- proach. In this paper, we will present a procedure to obtain the full real space interatomic force constant matrices using the direct supercell approach. It is based on the observa- tion that the planar force constants are in fact linear combinations of these matrix elements. Therefore from a limited set of planar force constants for some high- symmetry directions (which can be easily and accurately determined using supercells), the three-dimensional force constant matrices can be constructed. These force con- stants of the perfect crystal are particularly useful in the study of the dynamical properties of other mixed systems [2,3, 15]. In addition, phonon frequencies associated with any wave vector in the Brillouin zone can be obtained from these force constants. The method requires only the standard total energy codes with 10-20 atoms per cell and can be easily applied to both insulators and metals. One distinct advantage of this method over the linear response theory [8] is its ability to study the anharmonic effects as well. In the force calculation, the relation be- tween the force and the displacement may not be exactly linear, nor along the same direction. Although only the harmonic term will be discussed in the present Letter, in- formation on anharmonicity can be readily available if higher-order terms are kept in the planar force calcula- tion. Following similar procedures outlined below, infor- mation on the cubic or higher-order force constants (which will be third- or higher-rank tensors) can be ob- tained. This will be the subject of further investigations. In the harmonic approximation, the energy change re- sulting from small displacements of atoms is usually writ- ten in the general form Uh„m = — g u'(R). D'~(R — R') u~(R'), 1 2 R, R', a, p where u (R) is the deviation from equilibrium of atom a in the unit cell associated with lattice vector R, and D'~(R — R') is the force constant matrix connecting atom a in unit cell R and atom p in unit cell R'. (a and p are indices of atoms in the basis. ) There are some gen- eral symmetries that must be obeyed by the matrices D'~(R) [16], yet it has not been computationally feasible to calculate these matrices by directly evaluating the LDA total energies for a series of geometries with isolat- 1992 The American Physical Society 2799