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