PHYSICAL REVIE%' B VOLUME 41, NUMBER 5 15 FEBRUARY 1990-I Fast iterative diagonalixation of nonlocal pseudopotential Hamiltonians using the fast Fourier transformation S. Goedecker and K. Maschke Institut de Physique Appliquee, Departement de Physique, Ecole Polytechnique Federale de Lausanne, PHB-Ecublens, CH-1015 Lausanne, Switzerland (Received 18 September 1989) Many electronic-structure calculations are using nonlocal pseudopotentials in connection with iterative diagonalization methods. The most time-consuming steps in those calculations are the cal- culations of the n' matrix elements and the subsequent diagonalization of this matrix, which re- quires the repeated calculation of matrix vector products. Using conventional matrix multiplication this takes an order of n2 operations. We present a new technique which significantly reduces the computational effort for k points of high symmetry. Since the Hamiltonian is defined as an opera- tor, no individual matrix elements are required. Furthermore, it allows the use of the fast Fourier transform to do the matrix times vector multiplication. INTRODUCTION It is well known that the matrix multiplication Hx, where H is the n X n Hamiltonian matrix expressed in plane waves, can be done with the help of a fast Fourier transformation with n 1og2n operations if the potential is local. ' This method takes advantage of the fact that the kinetic energy is diagonal in Fourier space whereas the potential energy is diagonal in real space. To calculate the potential energy, one therefore first transforms into real space, multiplies with the potential, and then trans- forrns back. For large n, this is significantly faster than conventional matrix multiplication. Unfortunately, no exact method has hitherto been known which allows the use of the fast Fourier transformation in connection with nonlocal potentials. This has prevented the use of nonse- parable nonlocal potentials in large-scale computations such as molecular-dynamics simulations. THK ALGORITHM FOR NONLOCAL POTENTIALS The nonlocal part of the pseudopotential Hamiltonian such as proposed by Bachelet, Hamann, and Schluter is of the form calculation will take an order of n log2n operations. We now show how the decomposition of the wave function 4 into the angular momentum components can be eSciently implemented. THK DECOMPOSITION OF THK WAVE FUNCTIONS INTO ANGULAR MOMENTUM COMPONENTS We start with the well-known expansion QO I e'"'= g i'j((kr) g F; (8„$, )I'( (8",P"), I=O m= — I where k is the sum of the Bloch wave vector and a reciprocal-lattice vector. Let us furthermore expand the spherical Bessel functions into an arbitrary set of func- tions f(„„(r }: A, , K We have chosen to index f by two indices for reasons which will become apparent later on. Inserting (2) into (I}we obtain 1=0 m= — l A, , K max X Y(' (8„$, } . (3) 1=0 where I, „ is typically 2. The basic idea is now the fol- lowing. To calculate HV, where 4= g(, c "e'"'", one first splits up 4 into components corresponding to the different angular moments and a remainder: 4 I — p+4 ) — ]+ ' +4 I — I +4, . Then one calcu- max lates 01 — 0 41 — & ... 41 — I 4, in real space and multi- max plies each part with the corresponding loca1 potential V((r). Finally one transforms back into Fourier space. Since the transformation from real space into Fourier space can be done with n logan operations, this part of the @I(r) = X &(. .. (, A, , K, m where (4) P(. . (, k We can now interpret the coefficients o;(, „, 1'(' (8", P") as the coeScients of a matrix whose columns are labeled by k and whose rows are labeled by A„(r, l, m. This matrix does the transformation from the basis set e'"' to the basis set f(„„(r)I'( (8„$, ). For arbitrary 4 we can therefore write 41 3230 1990 The American Physical Society