1!--------------------------------------------------------------------
2subroutine dvex(nu,dvy)
3   !--------------------------------------------------------------------
4   !
5   USE kinds,     ONLY: DP
6   USE constants, ONLY: electron_charge_squared => e2
7   USE ld1inc,    ONLY: nwf, &
8                        psi, &
9                        get_angular_momentum => ll, &
10                        occupation => oc,&
11                        sl3, &
12                        nspin, &
13                        get_spin => isw, &
14                        grid
15   use radial_grids, only: ndmx, hartree
16
17   implicit none
18   !
19   ! I/O variables
20   !
21   integer,intent(in) :: nu
22   real (DP),intent(out) :: dvy(ndmx)
23
24
25
26   !
27   ! local variables
28   !
29   integer :: i, mu, l0, l1, l2, l3
30   real (DP) :: wrk(ndmx), wrk1(ndmx)
31   real (DP) :: fac, sss, &
32              ocs, &
33              doc, &
34              half
35   !
36   do i=1,grid%mesh
37      dvy(i)=0.0d0
38   end do
39   l1 = get_angular_momentum(nu)
40   half = 2.d0 * l1 + 1.d0
41
42
43   ! iterate over all the KS wavefunctions
44   ! get_spin(i) gives the spin of the i wavefunction
45   do mu = 1,nwf
46      !
47      ! only wfc with the same spin contribute to exchange term
48      !
49      ! check for the spin of the wavefunction
50      if (get_spin(mu) /= get_spin(nu) ) cycle  ! skip to next
51      ocs = occupation(mu) * (0.5d0 * nspin)
52!      write (*,*) mu, oc(mu), ocs
53      if ( mu == nu ) then
54         doc = 0.d0
55         !if(AND((l1 /= 0), (ocs > 0.d0))) then
56         if((l1 /= 0).AND.(ocs > 0.d0)) then
57           i = int(ocs)
58           doc = (i*(2.d0*ocs-i-1.d0)/(half-1.d0) - ocs*ocs/half) * half/ocs
59         end if
60         ocs = ocs + doc
61      !   if (doc /= 0.d0) write (*,*) "DOC ",nu, doc
62      end if
63!
64      l2 = get_angular_momentum(mu)
65      l0 = abs(l1 - l2)
66      wrk = psi(:,1,mu)*psi(:,1,nu)
67      ! do i = 1,grid%mesh
68
69      !    wrk(i) = psi(i,1,mu)*psi(i,1,nu) ! \varphi^{*}_{mu \sigma}\varphi^{*}_{nu \sigma}
70      ! end do
71      do l3 = l0, l1 + l2
72         sss = sl3(l1,l2,l3) ! some kind of symbol
73!         write (*,*) l1,l2,l3,sss
74         if (abs(sss).gt.1.0d-10) then
75            ! call the hartree utility function
76            ! from radial_grids
77            ! the output goes to vh
78            call hartree(l3,l1+l2+2,grid%mesh, grid, wrk, vh = wrk1)
79
80            fac =-electron_charge_squared*ocs*sss/2.0d0
81
82            do i=1,grid%mesh
83               dvy(i)= dvy(i) + fac*wrk1(i)*psi(i,1,mu)
84
85            end do
86         end if
87      end do
88
89
90
91
92!- spurious hartree part
93
94      if (mu == nu ) then
95          call hartree(0,2,grid%mesh,grid,wrk,wrk1)
96          fac = doc  * electron_charge_squared
97          do i=1,grid%mesh
98             dvy(i)= dvy(i) + fac*wrk1(i)*psi(i,1,mu)
99          end do
100      end if
101
102   end do
103
104end subroutine dvex
105