1subroutine fegrad(enrgy,bantot,rprimd,ngkpt,kpt,nkpt,shiftk,ikndx,indxkbnd,ncband,egrad)
2implicit none
3integer :: bantot,ngkpt(3),nkpt,ncband
4integer :: indxkbnd(nkpt),ikndx(ngkpt(1),ngkpt(2),ngkpt(3))
5double precision :: twopi
6double precision :: enrgy(bantot),rprimd(3,3),shiftk(3),kpt(3,nkpt)
7double precision :: egrad(3,bantot)
8double precision :: engrid(ngkpt(1),ngkpt(2),ngkpt(3))
9double precision :: fx(ngkpt(1),ngkpt(2),ngkpt(3))
10double precision :: fy(ngkpt(1),ngkpt(2),ngkpt(3))
11double precision :: fz(ngkpt(1),ngkpt(2),ngkpt(3))
12double precision :: fxy(ngkpt(1),ngkpt(2),ngkpt(3))
13double precision :: fxz(ngkpt(1),ngkpt(2),ngkpt(3))
14double precision :: fyz(ngkpt(1),ngkpt(2),ngkpt(3))
15double precision :: fxyz(ngkpt(1),ngkpt(2),ngkpt(3))
16double precision :: gridx(ngkpt(1)),gridy(ngkpt(2)),gridz(ngkpt(3))
17double precision :: xx,yy,zz,grk(3)
18integer :: ix,iy,iz,ii,jj,ie,ikk(3),iband,ikpt
19
20! Unit cartesian vectors x(ii,1), x(ii,2), x(ii,3)
21! Reciprocal lattice vector jj G(ii,jj)
22! Direct lattice vector jj R(ii,jj)
23! Wave vector expressed in either cartesian coordiantes or reciprocal lattice coordinates
24! k(ii) = x(ii,1) p(1) + x(ii,2) p(2) + x(ii,3) p(3)
25!       = G(ii,1) q(1) + G(ii,2) q(2) + G(ii,3) q(3)
26! sum_ii R(ii,jj)*G(ii,kk) = 2 pi delta(jj,kk)
27! implies q(jj) = sum_ii R(ii,jj)*k(ii)/(2 pi)
28! which implies (partial/partial p(ii)) q(jj) = sum_kk x(kk,ii)*R(kk,jj)/(2 pi)
29!                                             = R(ii,jj)/(2 pi)
30! Energy known on q grid E(q(1),q(2),q(3))
31! Want cartesian gradient of E
32! grad E(ii) = (partial/partial p(ii)) E(q(1),q(2),q(3))
33! = sum_jj [(partial/partial p(ii)) q(jj)] (partial/partial q(jj)) E(q(1),q(2),q(3))
34! = sum_jj R(ii,jj) (partial/partial q(jj)) E(q(1),q(2),q(3)) / (2 pi)
35
36twopi=2.d0*acos(-1.d0)
37do ix=1,ngkpt(1)
38  xx=dble(ix)/dble(ngkpt(1))+shiftk(1)/dble(ngkpt(1))
39  xx=mod(xx,1.d0)
40  if (xx.lt.0.d0) xx=xx+1.d0
41  xx=1.d0-mod(1.d0-xx,1.d0)
42  xx=xx-0.5d0
43  gridx(ix)=xx
44enddo
45do iy=1,ngkpt(2)
46  yy=dble(iy)/dble(ngkpt(2))+shiftk(2)/dble(ngkpt(2))
47  yy=mod(yy,1.d0)
48  if (yy.lt.0.d0) yy=yy+1.d0
49  yy=1.d0-mod(1.d0-yy,1.d0)
50  yy=yy-0.5d0
51  gridy(iy)=yy
52enddo
53do iz=1,ngkpt(3)
54  zz=dble(iz)/dble(ngkpt(3))+shiftk(3)/dble(ngkpt(3))
55  zz=mod(zz,1.d0)
56  if (zz.lt.0.d0) zz=zz+1.d0
57  zz=1.d0-mod(1.d0-zz,1.d0)
58  zz=zz-0.5d0
59  gridz(iz)=zz
60enddo
61do iband=1,ncband
62  do ix=1,ngkpt(1)
63  do iy=1,ngkpt(2)
64  do iz=1,ngkpt(3)
65    ikpt=ikndx(ix,iy,iz)
66    engrid(ix,iy,iz)=enrgy(indxkbnd(ikpt)+iband)
67  enddo
68  enddo
69  enddo
70  call xyzspline(engrid,gridx,gridy,gridz,ngkpt(1),ngkpt(2),ngkpt(3),1.d0,1.d0,1.d0,fx,fy,fz,fxy,fxz,fyz,fxyz)
71  do ikpt=1,nkpt
72    do ii=1,3
73      ikk(ii)=nint((kpt(ii,ikpt)+0.5d0)*dble(ngkpt(ii))-shiftk(ii))
74    enddo
75    grk=(/fx(ikk(1),ikk(2),ikk(3)),fy(ikk(1),ikk(2),ikk(3)),fz(ikk(1),ikk(2),ikk(3))/)
76    ie=indxkbnd(ikpt)+iband
77    do ii=1,3
78      egrad(ii,ie)=0.d0
79      do jj=1,3
80        egrad(ii,ie)=egrad(ii,ie)+rprimd(jj,ii)*grk(jj)/twopi
81      enddo
82    enddo
83  enddo
84enddo
85
86return
87end
88