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