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