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