1!----------------------------------------------------------------------- 2 subroutine kin_e_density & 3 (ndm, mesh, nwf, ll, oc, chi, r, r2, dx, tau ) 4!----------------------------------------------------------------------- 5! 6! calculate the (spherical) kinetic energy density 7! tau = sum_nwf 0.5*|\nabla phi|^2 in Cartesian coordinate 8! = 0.5*sum_nwf occ(nwf)/4pi*{[d(chi/r)/dr]^2 + (chi/r)^2*ll*(ll+1)/r^2} 9! = 0.5*sum_nwf occ(nwf)/4pi*[(1/r*dchi-chi/r^2)^2+(chi/r)^2*ll(ll+1)/r^2)] 10! in spherical coordinate 11! temporary - must be rewritten in a cleaner way 12! 13 implicit none 14! input 15 integer ndm, mesh, nwf,ll(nwf) 16 real*8 oc(nwf), chi(ndm,2,nwf),r2(ndm),r(ndm),dx 17! output 18 real*8 tau(ndm,2) 19! local 20 logical tspin 21 integer i,n 22 real*8 z,full,half,rho_tot, fourpi, & 23 & ocup,ocdn,temp,nll, corr 24 parameter ( fourpi = 4.d0*3.141592653589793d+00 ) 25 real*8 dchi(ndm) 26 tspin = .false. 27! 28 tau (:,:)=0.0d0 29! 30 do n=1,nwf 31 half=2*ll(n)+1 32 full=half+half 33 nll=ll(n)*(ll(n)+1) 34 call deriv5pt(mesh,dx,r,chi(1,1,n),dchi) 35 if( oc(n).le.half) then 36 ocup=oc(n) 37 ocdn=0.d0 38 else 39 tspin = .true. 40 ocup=half 41 ocdn=oc(n)-half 42 end if 43! 44 do i=1,mesh 45! kinetic energy density 46 temp=(dchi(i)-chi(i,1,n)/r(i))**2.0d0 & 47 & + chi(i,1,n)**2.0d0/r2(i)*nll 48 tau(i,1)=tau(i,1)+ocup*temp 49 tau(i,2)=tau(i,2)+ocdn*temp 50 end do 51 end do 52! 53 do i=1,mesh 54 tau(i,:)=tau(i,:)/fourpi*0.5d0/r2(i) 55 end do 56 return 57 end 58! 59 subroutine deriv5pt(mesh,dx,r,v,dv) 60! 61! numerical derivative using 5-point formula - error: O(dx^5) 62! Assumes a logarithmic grid r_i = r_0 exp((i-1)dx) 63! 64 implicit none 65! input 66 integer mesh 67 real*8 v(mesh), r(mesh), dx 68! output: dv = dv/dr = (1/r) dv/dx 69 real*8 dv(mesh) 70! 71 integer i 72! 73! 74 dv(1)=(-25.d0*v(1)+48.d0*v(2)-36.d0*v(3)+16.d0*v(4) & 75 & -3.d0*v(5))/(12.d0*dx*r(1)) 76 dv(2)=(-3.d0*v(1)-10.d0*v(2)+18.d0*v(3)- 6.d0*v(4) & 77 & + v(5))/(12.d0*dx*r(2)) 78! 79 do i=3,mesh-2 80 dv(i)=(v(i-2)-8.d0*v(i-1)+8.d0*v(i+1)-v(i+2))/(12.d0*dx*r(i)) 81 end do 82! 83 dv(mesh-1)=( 3.d0*v(mesh)+10.d0*v(mesh-1)-18.d0*v(mesh-2)+ & 84 & 6.d0*v(mesh-3)- v(mesh-4))/(12.d0*dx*r(mesh-1)) 85 dv(mesh)=( 25.d0*v(mesh) -48.d0*v(mesh-1)+36.d0*v(mesh-2)- & 86 & 16.d0*v(mesh-3)+3.d0*v(mesh-4))/(12.d0*dx*r(mesh)) 87! 88 return 89 end 90