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