1!--------------------------------------------------------------------- 2subroutine dfx_new(dchi0, vx) 3!--------------------------------------------------------------------- 4! OEP calculation of the exchange potential 5! -------------------------------------------------------------------- 6#undef DEBUG 7 8 use constants, only: pi 9 use kinds, only: DP 10 use ld1_parameters, only : nwfx 11 use ld1inc, only: nwf, nspin, oc, rho, psi, isw, grid 12 use radial_grids, only: ndmx, hartree 13 implicit none 14 ! 15 ! I/O variables 16 ! 17 real(DP) :: dchi0(ndmx,nwfx) 18 real(DP),intent(out) :: vx(ndmx,2) ! the oep exchange only potential 19 ! for the two spin configurations 20 21 ! 22 ! local variables 23 ! 24 integer, parameter :: niterx = 12 ! 6, 12, 24 25 real(DP) :: drho0(ndmx,2),appchim1(ndmx,2), vslater(ndmx,2) 26 real(DP) :: int_0_inf_dr 27 real(DP) :: vvx(ndmx,2,niterx),drhox(ndmx,2,niterx), dvh(ndmx,2,niterx), & 28 dvh0(ndmx,2), aux(ndmx), drho1(ndmx,2), dvh1(ndmx,2) 29 real(DP) :: a(niterx,niterx), & 30 inva(niterx,niterx), & ! inverse of a, i guess 31 b1(niterx), b2(niterx), c, c1, work(niterx), x(niterx), uno 32 integer :: iwork(niterx), info, iterx 33 integer :: i, j, jter, k, nu, is 34 real(DP) :: third, fac, capel 35 logical :: first = .true. 36 save first 37 38!-set the left hand side of the equation 39 call drho0ofvx(drho0,dchi0) 40 41 do is=1,nspin 42 call hartree(0,2,grid%mesh,grid,drho0(1,is),dvh0(1,is)) 43 end do 44 45 aux(1:grid%mesh) = drho0(1:grid%mesh,1)*dvh0(1:grid%mesh,1) 46 if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + & 47 drho0(1:grid%mesh,2)*dvh0(1:grid%mesh,2) ) 48 c = int_0_inf_dr(aux,grid,grid%mesh,2) 49 50 if (first) then 51! first = .false. 52!-slater exchange potential 53 do is=1,nspin 54 vslater(1:grid%mesh,is) = 0.d0 55 do nu=1,nwf 56 if (isw(nu) == is) vslater(1:grid%mesh,is) = vslater(1:grid%mesh,is) + & 57 oc(nu)*psi(1:grid%mesh,1,nu)*dchi0(1:grid%mesh,nu) 58 end do 59 vslater(1:grid%mesh,is) = vslater(1:grid%mesh,is) / rho(1:grid%mesh,is) 60! write (*,*) vslater(1:grid%mesh,is) 61 end do 62 else 63 vslater(1:grid%mesh,1:nspin) = vx(1:grid%mesh,1:nspin) 64 end if 65!- is a reasonable starting guess 66 call drhoofv(drho1,vslater) 67 68 drho1(1:grid%mesh,1:nspin) = drho1(1:grid%mesh,1:nspin) - drho0(1:grid%mesh,1:nspin) 69 do is=1,nspin 70 call hartree(0,2,grid%mesh,grid,drho1(1,is),dvh1(1,is)) 71 end do 72 73 aux(1:grid%mesh) = drho1(1:grid%mesh,1) * dvh1(1:grid%mesh,1) 74 if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + & 75 drho1(1:grid%mesh,2) * dvh1(1:grid%mesh,2) ) 76 c1 = int_0_inf_dr(aux,grid,grid%mesh,2) 77 78#ifdef DEBUG 79 write (*,'(a,2f16.10)') "C, C1 ", c, c1 80#endif 81 82!- simple Thomas-Fermi approximation to \chi^-1 83 third = 1.0d0/3.0d0 84 fac = -(0.75d0*pi)**(2.0d0/3.0d0) 85 86 do is =1, nspin 87 appchim1(1:grid%mesh,is) = fac/(grid%r(1:grid%mesh)* & 88 (nspin*rho(1:grid%mesh,is)*grid%r(1:grid%mesh))**third) 89 end do 90 91 drhox(1:grid%mesh,1:nspin,1) = drho1(1:grid%mesh,1:nspin) 92 93 if (c1 < 1.d-12 ) then 94 vx(1:grid%mesh,1:nspin) = vslater(1:grid%mesh,1:nspin) 95 return 96 end if 97!- ITERATE ! 98 do iterx =1,niterx 99!- set a new normalized correction vector vvx = chim1*drho/norm 100 101 vvx(1:grid%mesh,1:nspin,iterx) = appchim1(1:grid%mesh,1:nspin) * & 102 drhox(1:grid%mesh,1:nspin,iterx) 103 do is=1,nspin 104 call hartree(0,2,grid%mesh,grid,vvx(1,is,iterx),dvh(1,is,iterx)) 105 end do 106 107 aux(1:grid%mesh) =vvx(1:grid%mesh,1,iterx) * dvh(1:grid%mesh,1,iterx) 108 if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + & 109 vvx(1:grid%mesh,2,iterx) * dvh(1:grid%mesh,2,iterx) ) 110 capel = int_0_inf_dr(aux,grid,grid%mesh,2) 111#ifdef DEBUG 112 write (*,*) "norm ", capel 113#endif 114 if (capel > 0) then 115 capel = 1.d0/sqrt(capel) 116 vvx(1:grid%mesh,1:nspin,iterx) = vvx(1:grid%mesh,1:nspin,iterx) * capel 117 end if 118!- compute the corresponding drho 119 call drhoofv(drhox(1,1,iterx),vvx(1,1,iterx) ) 120 do is =1,nspin 121 call hartree(0,2,grid%mesh,grid,drhox(1,is,iterx),dvh(1,is,iterx)) 122 end do 123 124 aux(1:grid%mesh) = drhox(1:grid%mesh,1,iterx) * dvh1(1:grid%mesh,1) 125 if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + & 126 drhox(1:grid%mesh,2,iterx) * dvh1(1:grid%mesh,2) ) 127 b1(iterx) = int_0_inf_dr(aux,grid,grid%mesh,2) 128 129 aux(1:grid%mesh) = drho1(1:grid%mesh,1) * dvh(1:grid%mesh,1,iterx) 130 if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + & 131 drho1(1:grid%mesh,2) * dvh(1:grid%mesh,2,iterx) ) 132 b2(iterx) = int_0_inf_dr(aux,grid,grid%mesh,2) 133 134 do jter =1,iterx 135 136 aux(1:grid%mesh) = drhox(1:grid%mesh,1,iterx) * dvh(1:grid%mesh,1,jter) 137 if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + & 138 drhox(1:grid%mesh,2,iterx)*dvh(1:grid%mesh,2,jter) ) 139 140 a(iterx,jter) = int_0_inf_dr(aux,grid,grid%mesh,2) 141 142 aux(1:grid%mesh) = drhox(1:grid%mesh,1,jter) * dvh(1:grid%mesh,1,iterx) 143 if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + & 144 drhox(1:grid%mesh,2,jter)*dvh(1:grid%mesh,2,iterx) ) 145 a(jter,iterx) = int_0_inf_dr(aux,grid,grid%mesh,2) 146 147 end do 148 149 capel = 0.d0 150 151 do i=1,iterx 152 do j=1,iterx 153 capel = capel + abs(a(i,j)-a(j,i)) 154 a(i,j) = 0.5d0*(a(i,j)+a(j,i)) 155 a(j,i) = a(i,j) 156 end do 157 end do 158 159 160#ifdef DEBUG 161 write (*,'(a,12f16.10)') "CAPEL a ", capel 162#endif 163 164 inva = a 165 166 167 168 CALL DSYTRF('U', iterx, inva, niterx, iwork, work, niterx,info) 169 if (info.ne.0) stop 'factorization' 170 CALL DSYTRI( 'U', iterx, inva, niterx, iwork, work, info ) 171 if (info.ne.0) stop 'DSYTRI' 172 forall ( i=1:iterx, j=1:iterx, j>i ) inva(j,i) = inva(i,j) 173 174! write (*,*) "INVA " 175! write (*,'(12f16.10)') inva 176 capel = 0.d0 177 do i=1,iterx 178 do j=1,iterx 179 uno = 0.d0 180 do k=1,iterx 181 uno = uno + a(i,k)*inva(k,j) 182 end do 183 if (i.eq.j) uno = uno - 1.d0 184 capel = capel + abs(uno) 185 end do 186 end do 187#ifdef DEBUG 188 write (*,'(a,12f16.10)') "CAPEL uno", capel 189#endif 190 191 x = 0.d0 192 capel = c1 193 do i=1,iterx 194 do j=1,iterx 195 x(i) = x(i) - inva(i,j) * 0.5d0*(b1(j)+b2(j)) 196 end do 197 capel = capel + x(i) * (b1(i)+b2(i)) 198 do j =1,i 199 capel = capel + x(i)*a(i,j)*x(j) 200 end do 201 do j =1,i-1 202 capel = capel + x(i)*a(j,i)*x(j) 203 end do 204 end do 205! write (*,'(a,12f16.10)') "X ", x 206#ifdef DEBUG 207 write (*,*) "capel ", capel 208#endif 209 210 211 vx(1:grid%mesh,1:nspin) = vslater(1:grid%mesh,1:nspin) 212 do j=1,iterx 213 vx(1:grid%mesh,1:nspin) = vx(1:grid%mesh,1:nspin) + x(j) * vvx(1:grid%mesh,1:nspin,j) 214 end do 215 216 if (iterx.eq.niterx) return 217 218 call drhoofv(drhox(1,1,iterx+1), vx ) 219 drhox(1:grid%mesh,1:nspin,iterx+1) = drhox(1:grid%mesh,1:nspin,iterx+1) - & 220 drho0(1:grid%mesh,1:nspin) 221 do is=1,nspin 222 call hartree(0,2,grid%mesh,grid,drhox(1,1,iterx+1),dvh(1,1,iterx+1)) 223 end do 224 225 aux(1:grid%mesh) = drhox(1:grid%mesh,1,iterx+1)*dvh(1:grid%mesh,1,iterx+1) 226 if (nspin==2) aux(1:grid%mesh) = 2.d0 * ( aux(1:grid%mesh) + & 227 drhox(1:grid%mesh,1,iterx+1)*dvh(1:grid%mesh,1,iterx+1)) 228 capel = int_0_inf_dr(aux,grid,grid%mesh,2) 229#ifdef DEBUG 230 write (*,*) "capel-check ", capel 231#endif 232! 233 end do 234! 235 return 236end subroutine dfx_new 237