1!
2!--------------------------------------------------------------------
3subroutine add_exchange ( energy )
4!--------------------------------------------------------------------
5#undef DEBUG
6   !
7   use io_global, only: stdout
8   use kinds,  only: DP
9   use constants, only: e2
10   use ld1_parameters, only : nwfx
11   use ld1inc, only: nwf, oc, psi, vx, sl3, ll, nn, enl, el, &
12                     isw, nspin, enzero, grid
13   use radial_grids, only: ndmx, hartree
14   implicit none
15   !
16   ! I/O variables
17   !
18   real (DP) :: energy
19   !
20   ! local variables
21   !
22   integer :: i, l0, l1, l2, l3, nu, mu, nst, is, half
23   real (DP) :: ex_hf, ocs, doc, sss, fac, sxc, sxc1
24   real (DP) :: wrk(ndmx), wrk1(ndmx), wrk2(ndmx), int_0_inf_dr, enzhf(nwfx)
25!
26   ex_hf = 0.0
27   do nu=1,nwf
28      is = isw(nu)
29      l1 = ll(nu)
30      half = 2.d0 * l1 + 1.d0
31      sxc = 0.0d0
32      do mu=1,nwf
33!
34! only wfc with the same spin contribute to exchange term
35!
36         if (isw(mu) /= is) cycle
37         ocs = oc(mu) * (0.5d0 * nspin )
38         if ( mu == nu ) then
39            doc = 0.d0
40            if( (l1 /= 0) .and. (ocs > 0.d0) ) then
41              i = int(ocs)
42              doc = (i*(2.d0*ocs-i-1.d0)/(half-1.d0) - ocs*ocs/half) * half/ocs
43            end if
44            ocs = ocs + doc
45!            if (doc /= 0.d0) write (*,*) "DOC ",nu, doc
46         end if
47
48         l2 = ll(mu)
49         l0=abs(l1-l2)
50         do i=1,grid%mesh
51            wrk(i) = psi(i,1,mu)*psi(i,1,nu)
52            wrk1(i)= 0.0d0
53         end do
54         do l3=l0,l1+l2
55            sss = sl3(l1,l2,l3)
56            if (abs(sss).gt.1.0d-10) then
57               call hartree(l3,l1+l2+2,grid%mesh,grid,wrk,wrk2)
58               fac = -e2*ocs*sss/2.0d0
59               do i=1,grid%mesh
60                  wrk1(i)= wrk1(i) + fac*wrk2(i)*wrk(i)
61               end do
62            end if
63         end do
64!- spurious hartree part
65          if (mu.eq.nu) then
66            call hartree(0,2,grid%mesh,grid,wrk,wrk2)
67            fac = doc*e2
68            do i=1,grid%mesh
69               wrk1(i) = wrk1(i) + fac*wrk2(i)*wrk(i)
70            end do
71          end if
72!
73         nst = 2 * min(l1,l2) + 2
74         sxc = sxc + int_0_inf_dr(wrk1, grid, grid%mesh, nst)
75      end do
76!
77      do i=1,grid%mesh
78         wrk1(i) = vx(i,is)*psi(i,1,nu)*psi(i,1,nu)
79      end do
80      sxc1 = int_0_inf_dr(wrk1,grid,grid%mesh,2*ll(nu)+2)
81      ex_hf = ex_hf + 0.5d0*oc(nu)*sxc
82      enzhf(nu)=sxc1-sxc
83      if(oc(nu)>0) enzero(is) = enzhf(nu)
84   end do
85
86   energy = energy + ex_hf
87
88#ifdef DEBUG
89   write (*,*) enzero
90   write(stdout,1100) (nn(nu),ll(nu),el(nu),oc(nu),enl(nu)-enzero(isw(nu)), &
91                  enl(nu)-enzhf(nu), nu=1,nwf)
92
931100 format(4x,2i2,5x,a2,'(',f5.2,')',2f11.4)
94#endif
95
96   return
97
98end subroutine add_exchange
99