1 2! Copyright (C) 2007 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl. 3! This file is distributed under the terms of the GNU General Public License. 4! See the file COPYING for license details. 5 6subroutine genvmatk(vmt,vir,ngp,igpig,wfmt,ld,wfgp,vmat) 7use modmain 8use modomp 9implicit none 10! arguments 11! the potential is multiplied by the radial integration weights in the 12! muffin-tin and by the characteristic function in the interstitial region 13real(8), intent(in) :: vmt(npcmtmax,natmtot),vir(ngtot) 14integer, intent(in) :: ngp(nspnfv),igpig(ngkmax,nspnfv) 15complex(8), intent(in) :: wfmt(npcmtmax,natmtot,nspinor,nstsv) 16integer, intent(in) :: ld 17complex(8), intent(in) :: wfgp(ld,nspinor,nstsv) 18complex(8), intent(out) :: vmat(nstsv,nstsv) 19! local variables 20integer ist,jst,ispn,jspn 21integer is,ias,npc,igp,nthd 22! automatic arrays 23complex(8) wfmt1(npcmtmax),z(ngkmax) 24! allocatable arrays 25complex(8), allocatable :: wfir(:) 26! external functions 27complex(8), external :: zdotc 28! zero the matrix elements 29vmat(:,:)=0.d0 30!-------------------------! 31! muffin-tin part ! 32!-------------------------! 33call holdthd(nstsv,nthd) 34!$OMP PARALLEL DO DEFAULT(SHARED) & 35!$OMP PRIVATE(wfmt1,ispn,ias) & 36!$OMP PRIVATE(is,npc,ist) & 37!$OMP NUM_THREADS(nthd) 38do jst=1,nstsv 39 do ispn=1,nspinor 40 do ias=1,natmtot 41 is=idxis(ias) 42 npc=npcmt(is) 43! apply potential to wavefunction 44 wfmt1(1:npc)=vmt(1:npc,ias)*wfmt(1:npc,ias,ispn,jst) 45! compute the inner products 46 do ist=1,jst 47 vmat(ist,jst)=vmat(ist,jst)+zdotc(npc,wfmt(:,ias,ispn,ist),1,wfmt1,1) 48 end do 49 end do 50 end do 51end do 52!$OMP END PARALLEL DO 53call freethd(nthd) 54!---------------------------! 55! interstitial part ! 56!---------------------------! 57call holdthd(nstsv,nthd) 58!$OMP PARALLEL DEFAULT(SHARED) & 59!$OMP PRIVATE(wfir,z,ispn,jspn) & 60!$OMP PRIVATE(igp,ist) & 61!$OMP NUM_THREADS(nthd) 62allocate(wfir(ngtot)) 63!$OMP DO 64do jst=1,nstsv 65 do ispn=1,nspinor 66 jspn=jspnfv(ispn) 67! Fourier transform wavefunction to real-space 68 wfir(:)=0.d0 69 do igp=1,ngp(jspn) 70 wfir(igfft(igpig(igp,jspn)))=wfgp(igp,ispn,jst) 71 end do 72 call zfftifc(3,ngridg,1,wfir) 73! apply potential to wavefunction 74 wfir(1:ngtot)=vir(1:ngtot)*wfir(1:ngtot) 75! Fourier transform to G+p-space 76 call zfftifc(3,ngridg,-1,wfir) 77 do igp=1,ngp(jspn) 78 z(igp)=wfir(igfft(igpig(igp,jspn))) 79 end do 80 do ist=1,jst 81! compute inner product 82 vmat(ist,jst)=vmat(ist,jst)+zdotc(ngp(jspn),wfgp(:,ispn,ist),1,z,1) 83 end do 84 end do 85end do 86!$OMP END DO 87deallocate(wfir) 88!$OMP END PARALLEL 89call freethd(nthd) 90! lower triangular part 91do ist=1,nstsv 92 do jst=1,ist-1 93 vmat(ist,jst)=conjg(vmat(jst,ist)) 94 end do 95end do 96end subroutine 97 98