1subroutine mkxclda(iband,ikpt,ngfft,vxc,vol,ihlf,nkpt,igndx,igmx,igmn,natom,exc) 2implicit none 3integer :: iband,ikpt,ngfft(3),ihlf(nkpt),igmn(3),igmx(3),natom,nkpt 4integer :: igndx(nkpt,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3)) 5double precision :: exc,vxc(ngfft(1)*ngfft(2)*ngfft(3)),vol 6integer :: ir,ix,iy,iz,iatom 7double precision :: rred(3),orbden 8double complex :: twfval,pwfval,wfval 9 10 exc=0.d0 11 ir=0 12 do iz=0,ngfft(3)-1 13 do iy=0,ngfft(2)-1 14 do ix=0,ngfft(1)-1 15 ir=ir+1 16 rred=dble((/ix,iy,iz/))/dble(ngfft) 17 call pwrwf(rred,iband,ikpt,twfval,ihlf,igndx,igmx,igmn) 18! do iatom=1,natom 19! call prjrwf(rred,iatom,iband,ikpt,pwfval) 20! enddo 21! wfval=twfval+pwfval 22 wfval=twfval 23 orbden=dble(wfval*conjg(wfval)) 24 exc=exc+orbden*vxc(ir) 25 enddo 26 enddo 27 enddo 28 exc=exc*vol/(ngfft(1)*ngfft(2)*ngfft(3)) 29 30return 31end subroutine mkxclda 32 33!****************************************************************************** 34 35subroutine pwrwf(rr,iband,ikpt,wfval,ihlf,igndx,igmx,igmn) 36! find plane wave expansion of wavefunction at position rr 37! rr in reduced coordinates 38use wfkvars 39use geometry 40implicit none 41integer :: iband,ikpt,ipw,ihlf(nkpt),igg(3),ipw2,igmx(3),igmn(3) 42integer :: igndx(nkpt,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3)) 43double complex :: phase,wfval,twopii 44double precision :: pi,rr(3) 45parameter (pi=3.1415926535897932384626433832795028841972d0) 46 47wfval=(0.d0,0.d0) 48twopii=(0.d0,2.d0)*pi 49do ipw=1,npwarr(ikpt) 50 phase=dot_product(rr,kpt(:,ikpt)+kg(:,indxkpw(ikpt)+ipw))*twopii 51 wfval=wfval+cg(indxkcg(ikpt)+(iband-1)*npwarr(ikpt)+ipw)*exp(phase) 52enddo 53wfval=wfval/sqrt(vol) 54 55return 56end subroutine pwrwf 57 58!****************************************************************************** 59 60subroutine prjrwf(rr,iatom,iband,ikpt,wfval) 61! find PAW correction to wavefunction at position rr in reduced coordinates 62use wfkvars 63use pawvars 64use geometry 65implicit none 66double precision :: rr(3),rat(3),rrel(3),phival,tphival,rrmag,linterp 67integer :: iband,ikpt,nn,ll,mm,iorb,ii,jj,iatom,itpaw 68double complex :: wfval,ylm,spharm 69external linterp,spharm 70 71 itpaw=typat(iatom) 72 rrel=(/0.d0,0.d0,0.d0/) 73 do ii=1,3 74 do jj=1,3 75 rrel(ii)=rrel(ii)+rprimd(ii,jj)*(rr(jj)-xred(jj,iatom)) 76 enddo 77 enddo 78 rrmag=sqrt(dot_product(rrel,rrel)) 79 wfval=0.d0 80 if (rrmag.gt.rphigrid(itpaw,gridsize(itpaw,iphigrid(itpaw)))) return 81 do iorb=1,nlmn(itpaw) 82 nn=iorbno(itpaw,iorb) 83 ll=llorb(itpaw,iorb) 84 mm=mmorb(itpaw,iorb) 85 ylm=spharm(ll,mm,rrel) 86 tphival=linterp(tphi(itpaw,nn,:),rphigrid(itpaw,:),gridsizemax,rrmag) 87 phival=linterp(phi(itpaw,nn,:),rphigrid(itpaw,:),gridsizemax,rrmag) 88 wfval=wfval+projwf(iatom,iorb,ikpt,iband)*ylm*(phival-tphival)/rrmag 89 enddo 90 91return 92end subroutine prjrwf 93 94!****************************************************************************** 95 96subroutine prjrtwf(rr,iatom,iband,ikpt,wfval) 97! find PAW correction to wavefunction at position rr in reduced coordinates 98use wfkvars 99use pawvars 100use geometry 101implicit none 102double precision :: rr(3),rat(3),rrel(3),phival,tphival,rrmag,linterp 103integer :: iband,ikpt,nn,ll,mm,iorb,ii,jj,iatom,itpaw 104double complex :: wfval,ylm,spharm 105external linterp,spharm 106 107 itpaw=typat(iatom) 108 rrel=(/0.d0,0.d0,0.d0/) 109 do ii=1,3 110 do jj=1,3 111 rrel(ii)=rrel(ii)+rprimd(ii,jj)*(rr(jj)-xred(jj,iatom)) 112 enddo 113 enddo 114 rrmag=sqrt(dot_product(rrel,rrel)) 115 wfval=0.d0 116 if (rrmag.gt.rphigrid(itpaw,gridsize(itpaw,iphigrid(itpaw)))) return 117 do iorb=1,nlmn(itpaw) 118 nn=iorbno(itpaw,iorb) 119 ll=llorb(itpaw,iorb) 120 mm=mmorb(itpaw,iorb) 121 ylm=spharm(ll,mm,rrel) 122 tphival=linterp(tphi(itpaw,nn,:),rphigrid(itpaw,:),gridsizemax,rrmag) 123 wfval=wfval+projwf(iatom,iorb,ikpt,iband)*ylm*tphival/rrmag 124 enddo 125 126return 127end subroutine prjrtwf 128 129!****************************************************************************** 130 131subroutine prjrpwf(rr,iatom,iband,ikpt,wfval) 132! find PAW correction to wavefunction at position rr in reduced coordinates 133use wfkvars 134use pawvars 135use geometry 136implicit none 137double precision :: rr(3),rat(3),rrel(3),phival,tphival,rrmag,linterp 138integer :: iband,ikpt,nn,ll,mm,iorb,ii,jj,iatom,itpaw 139double complex :: wfval,ylm,spharm 140external linterp,spharm 141 142 itpaw=typat(iatom) 143 rrel=(/0.d0,0.d0,0.d0/) 144 do ii=1,3 145 do jj=1,3 146 rrel(ii)=rrel(ii)+rprimd(ii,jj)*(rr(jj)-xred(jj,iatom)) 147 enddo 148 enddo 149 rrmag=sqrt(dot_product(rrel,rrel)) 150 wfval=0.d0 151 if (rrmag.gt.rphigrid(itpaw,gridsize(itpaw,iphigrid(itpaw)))) return 152 do iorb=1,nlmn(itpaw) 153 nn=iorbno(itpaw,iorb) 154 ll=llorb(itpaw,iorb) 155 mm=mmorb(itpaw,iorb) 156 ylm=spharm(ll,mm,rrel) 157 phival=linterp(phi(itpaw,nn,:),rphigrid(itpaw,:),gridsizemax,rrmag) 158 wfval=wfval+projwf(iatom,iorb,ikpt,iband)*ylm*phival/rrmag 159 enddo 160 161return 162end subroutine prjrpwf 163 164