1subroutine plus_u_setup(natih, lsr) 2! 3! Add additional +U orbitals (if DFT+U) to the full list of projectors 4! 5! 6 USE kinds, ONLY : DP 7 USE constants, ONLY : rytoev 8 use noncollin_module, ONLY : noncolin 9 USE ldaU, ONLY : lda_plus_U, lda_plus_u_kind, U_projection, & 10 Hubbard_lmax, Hubbard_l, Hubbard_U, Hubbard_alpha, & 11 Hubbard_J0, Hubbard_beta 12 use atom, ONLY : rgrid 13 USE scf, ONLY : rho 14 use radial_grids, ONLY : ndmx 15 USE ions_base, ONLY : nat, ityp, ntyp => nsp, atm 16 USE cell_base, ONLY : alat 17 use uspp_param, only : nhm, upf 18 USE io_global, ONLY : stdout 19 USE cond, ONLY : norbs, nocrosl, noinss, nocrosr, tblms, taunews, & 20 nenergy, earr, nrzs, zs, tran_tot, norbf, nbrx, & 21 cross, zpseus, zpseus_nc, betars, iofspin 22 implicit none 23 24 integer :: lsr, iorb, iorb1, it, iwfc, iwfc1, mesh, i, ipol, ldim, & 25 norbs_new, nocrosl_new, noinss_new, nocrosr_new, na, & 26 natih(2,norbs), lll, kkbeta 27 integer, allocatable :: ind(:,:), tblms_new(:,:), cross_new(:,:) 28 real(DP), parameter :: epswfc=1.d-4, eps=1.d-8 29 REAL(DP) :: r1, beta1, beta2, norm, ledge, redge 30 REAL(DP), ALLOCATABLE :: bphi(:,:), rsph(:), taunews_new(:,:), & 31 gi(:), zpseus_new(:,:,:) 32 33!-- 34! Some checks 35 36 if (.not.lda_plus_u) return 37 if (lda_plus_u_kind.eq.1) call errore('plus_u_setup','Full LDA+U not yet implemented',1) 38 39 WRITE( stdout, '(/,/,5x,"Simplified LDA+U calculation (l_max = ",i1, & 40 & ") with parameters (eV):")') Hubbard_lmax 41 WRITE( stdout, '(5x,A)') & 42 "atomic species L U alpha J0 beta" 43 DO it = 1, ntyp 44 IF ( Hubbard_U(it) /= 0.D0 .OR. Hubbard_alpha(it) /= 0.D0 .OR. & 45 Hubbard_J0(it) /= 0.D0 .OR. Hubbard_beta(it) /= 0.D0 ) THEN 46 WRITE( stdout,'(5x,a6,12x,i1,2x,4f9.4)') atm(it), Hubbard_L(it), & 47 Hubbard_U(it)*rytoev, Hubbard_alpha(it)*rytoev, & 48 Hubbard_J0(it)*rytoev, Hubbard_beta(it)*rytoev 49 END IF 50 END DO 51 52 if (U_projection.eq."pseudo") return 53 if (U_projection.ne."atomic") & 54 call errore('plus_u_setup','+U works only for U_projection=''pseudo'' or ''atomic'' ',1) 55 if (noncolin) call errore('plus_u_setup','+U for noncollinear case not yet implemented',1) 56 if (lsr.ne.2) call errore('plus_u_setup','+U atoms are allowed only in scatt. region',1) 57!-- 58 59 allocate ( gi(ndmx) ) 60 allocate ( bphi(nbrx,ntyp) ) 61 allocate ( rsph(ntyp) ) 62 allocate ( ind(2,norbs) ) 63 64 bphi(:,:) = 0.d0 65 rsph(:) = 0.d0 66 ind(:,:) = 0 67 68!-- 69! Calculate the total number of orbitals (beta + U WF's) 70! 71 72 noinss_new = noinss 73 norbs_new = norbs 74 75 iorb = 1 76 do while (iorb.le.norbs) 77 it = tblms(1,iorb) 78 if (Hubbard_U(it).ne.0.d0) then 79 iorb1 = iorb 80 do while (natih(1,iorb1).eq.natih(1,iorb)) 81 iorb1 = iorb1 + 1 82 enddo 83! 84! The last beta for the atom with U is provided with the 85! index of the 1st (iorb) and the last (iorb1) beta 86! 87 iorb1 = iorb1 - 1 88 ind(1,iorb1) = iorb 89 ind(2,iorb1) = iorb1 90!-- 91 ldim = 2*Hubbard_l(it)+1 92 noinss_new = noinss_new + ldim 93 norbs_new = norbs_new + ldim 94 iorb = iorb1 95 endif 96 iorb = iorb + 1 97 enddo 98!-- 99 100!-- 101! Determine the radii of atomic U WF's 102! 103 do it = 1, ntyp 104 if (Hubbard_U(it).ne.0.d0) then 105 do iwfc = 1, upf(it)%nwfc 106 if (upf(it)%lchi(iwfc).eq.Hubbard_l(it)) then 107 r1 = 0.d0 108 do i = 2, rgrid(it)%mesh 109 r1 = max(r1, ABS(upf(it)%chi(i,iwfc)/rgrid(it)%r(i))) 110 enddo 111 i = rgrid(it)%mesh 112 do while (abs(upf(it)%chi(i,iwfc)/rgrid(it)%r(i)).le.epswfc*r1) 113 i = i - 1 114 enddo 115 rsph(it) = rgrid(it)%r(i) / alat 116 mesh = i 117 endif 118 enddo 119 endif 120 enddo 121!-- 122 123!-- 124! Check that all +U orbitals are totally inside the scatt. region 125 126 i = 0 127 write(6,*) 'Scatt. region L = ', zs(nrzs+1) 128 do iorb = 1, norbs 129 if (ind(2,iorb).eq.iorb) then 130 it = tblms(1,iorb) 131 beta1 = taunews(3,iorb)-rsph(it) 132 beta2 = taunews(3,iorb)+rsph(it) 133 if (beta1.le.1.d-4.or.beta2.gt.zs(nrzs+1)-1.d-4) i = 1 134 endif 135 enddo 136 if (i.eq.1) call errore('plus_u_setup','some +U orbitals cross the boundary (not allowed) ...',1) 137!-- 138 139!-- 140! Calculate the integrals of betas with U atomic orbitals, 141! bphi(i) = \sum_j q_ij <beta_j|phi> 142! 143 do it = 1, ntyp 144 if (Hubbard_U(it).ne.0.d0) then 145 146 mesh = upf(it)%mesh 147 kkbeta = upf(it)%kkbeta 148 do iwfc = 1, upf(it)%nwfc 149 if (upf(it)%lchi(iwfc).eq.Hubbard_l(it)) then 150 do iorb = 1, upf(it)%nbeta 151 if (upf(it)%lll(iorb).eq.Hubbard_l(it)) then 152 gi(1:kkbeta)= upf(it)%beta(1:kkbeta,iorb) * & 153 upf(it)%chi (1:kkbeta,iwfc) 154 call simpson (kkbeta, gi, upf(it)%rab,bphi(iorb,it)) 155 endif 156 enddo 157 endif 158 enddo 159 gi(:) = 0.d0 160 do iorb = 1, upf(it)%nbeta 161 do iorb1 = 1, upf(it)%nbeta 162 gi(iorb) = gi(iorb) + upf(it)%qqq(iorb,iorb1)*bphi(iorb1,it) 163 enddo 164 enddo 165 bphi(1:upf(it)%nbeta,it) = gi(1:upf(it)%nbeta) 166 167 endif 168 enddo 169!-- 170 171!-- 172! Allocate the arrays with all the orbitals (beta + U WF's) 173! 174 allocate( taunews_new(4,norbs_new) ) 175 allocate( tblms_new(4,norbs_new) ) 176 allocate( cross_new(norbs_new, nrzs) ) 177 allocate( zpseus_new(2, norbs_new, norbs_new) ) 178 zpseus_new(:,:,:) = 0.d0 179!-- 180 181!-- 182! Set up new extended arrays (beta + U orbitals) 183! 184! iorb --> old list 185! iorb1 --> new list 186! 187! old list new list 188! 189! (1st atom beta ) (1st atom beta ) 190! iorb --> (2nd atom beta ) ( +U orbitals ) 191! ( ... ) iorb1 --> (2nd atom beta ) 192! ( +U orbitals ) 193! ( ... ) 194 iorb1 = 0 195 do iorb = 1, norbs 196 iorb1 = iorb1 + 1 197 it = tblms(1,iorb) 198 na = natih(1,iorb) 199!-- 200! setting up some beta arrays from old ones (just shifting) 201 202 do i = 1, 4 203 taunews_new(i,iorb1) = taunews(i,iorb) 204 enddo 205 do i = 1, 4 206 tblms_new(i,iorb1) = tblms(i,iorb) 207 enddo 208 do i = 1, nrzs 209 cross_new(iorb1,i) = cross(iorb,i) 210 enddo 211!-- 212 213!-- 214! beta-beta block of zpseu (again just shifting) 215! 216 do i = 1, norbs 217 if(natih(1,i).eq.na) then 218 zpseus_new(:,i-iorb+iorb1,iorb1) = zpseus(:,i,iorb) 219 endif 220 enddo 221!-- 222 223! entering into +U orbitals part (if any) 224 225 if (ind(2,iorb).eq.iorb) then 226 227 lll = Hubbard_l(it) 228 ldim = 2*lll + 1 229!-- 230! beta-beta additional block of zpseu 231! 232 do iwfc = ind(1,iorb), iorb 233 if (tblms(3,iwfc).eq.lll) then 234 do iwfc1 = ind(1,iorb), iorb 235 if (tblms(3,iwfc1).eq.lll) then 236 r1 = -2.d0*rho%ns(tblms(4,iwfc),tblms(4,iwfc1),iofspin,na) 237 if (tblms(4,iwfc).eq.tblms(4,iwfc1)) r1 = r1 + 1.d0 238 zpseus_new(1,iwfc-iorb+iorb1,iwfc1-iorb+iorb1) = & 239 zpseus_new(1,iwfc-iorb+iorb1,iwfc1-iorb+iorb1) + & 240 0.5d0*Hubbard_U(it)*bphi(tblms(2,iwfc),it)*bphi(tblms(2,iwfc1),it)*r1 241 endif 242 enddo 243 endif 244 enddo 245!-- 246 247!-- 248! beta-atomic block of zpseu 249! 250 lll = Hubbard_l(it) 251 do iwfc = ind(1,iorb), iorb 252 if (tblms(3,iwfc).eq.lll) then 253 do iwfc1 = 1, ldim 254 r1 = -2.d0*rho%ns(tblms(4,iwfc),iwfc1,iofspin,na) 255 if (tblms(4,iwfc).eq.iwfc1) r1 = r1 + 1.d0 256 zpseus_new(1,iwfc-iorb+iorb1,iorb1+iwfc1) = & 257 0.5d0*Hubbard_U(it)*bphi(tblms(2,iwfc),it)*r1 258 enddo 259 endif 260 enddo 261!-- 262 263!-- 264! atomic-beta block of zpseu 265! 266 do iwfc1 = 1, ldim 267 do iwfc = ind(1,iorb), iorb 268 zpseus_new(1,iorb1+iwfc1,iwfc-iorb+iorb1) = & 269 zpseus_new(1,iwfc-iorb+iorb1,iorb1+iwfc1) 270 enddo 271 enddo 272!-- 273 274!-- 275! atomic-atomic block of zpseu 276! 277 do iwfc = 1, ldim 278 do iwfc1 = 1, ldim 279 zpseus_new(1,iorb1+iwfc,iorb1+iwfc1) = & 280 - Hubbard_U(it) * rho%ns(iwfc,iwfc1,iofspin,na) 281 enddo 282 zpseus_new(1,iorb1+iwfc,iorb1+iwfc) = & 283 zpseus_new(1,iorb1+iwfc,iorb1+iwfc) + 0.5d0*Hubbard_U(it) 284 enddo 285!-- 286 287!-- 288! setting up some +U orbitals arrays from those of beta 289! 290 291 do iwfc = 1, ldim 292 293 iorb1 = iorb1 + 1 294 295 do i = 1, 3 296 taunews_new(i,iorb1) = taunews(i,iorb) 297 enddo 298 taunews_new(4,iorb1) = rsph(it) 299 300 tblms_new(1,iorb1) = tblms(1,iorb) 301 tblms_new(2,iorb1) = tblms(2,iorb) + 1 302 tblms_new(3,iorb1) = Hubbard_l(it) 303 tblms_new(4,iorb1) = iwfc 304 305 ledge = taunews(3,iorb) - rsph(it) 306 redge = taunews(3,iorb) + rsph(it) 307 do i = 1, nrzs 308 if (ledge.gt.zs(i+1).or.redge.lt.zs(i)) then 309 cross_new(iorb1,i)=0 310 else 311 cross_new(iorb1,i)=1 312 endif 313 enddo 314 315 enddo 316!-- 317 318 endif 319 320 enddo 321!-- 322 323!-- 324! Add the atomic radial WF's with U to the list betars 325! 326 do it = 1, ntyp 327 if (Hubbard_U(it).ne.0.d0) then 328 do iwfc = 1, upf(it)%nwfc 329 if (upf(it)%lchi(iwfc).eq.Hubbard_l(it)) then 330 betars(1:rgrid(it)%mesh,upf(it)%nbeta+1,it) = & 331 upf(it)%chi(1:rgrid(it)%mesh,iwfc) 332 endif 333 enddo 334 endif 335 enddo 336!-- 337 338!-- 339! Reallocate the orbital arrays with new dimensions and date 340! 341 deallocate( taunews ) 342 deallocate( tblms ) 343 deallocate( cross ) 344 deallocate( zpseus ) 345 346 noinss = noinss_new 347 norbs = norbs_new 348 norbf = norbs 349 350 allocate( taunews(4,norbs) ) 351 allocate( tblms(4,norbs) ) 352 allocate( cross(norbs, nrzs) ) 353 354 taunews(:,:) = taunews_new(:,:) 355 tblms(:,:) = tblms_new(:,:) 356 cross(:,:) = cross_new(:,:) 357 358 allocate( zpseus(2, norbs, norbs) ) 359 zpseus(:,:,:) = zpseus_new(:,:,:) 360!-- 361 deallocate( gi ) 362 deallocate( bphi ) 363 deallocate( rsph ) 364 deallocate( ind ) 365 366 deallocate( taunews_new ) 367 deallocate( tblms_new ) 368 deallocate( cross_new ) 369 deallocate( zpseus_new ) 370 371 return 372end subroutine plus_u_setup 373 374