TABLE OF CONTENTS
ABINIT/pawpuxinit [ Functions ]
NAME
pawpuxinit
FUNCTION
Initialize some starting values of several arrays used in PAW+U/+DMFT or local exact-exchange calculations A-define useful indices for LDA+U/local exact-exchange B-Compute overlap between atomic wavefunction C-Compute matrix elements of coulomb interaction (see PRB vol.52 5467) (angular part computed from Gaunt coefficients)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (BA,FJ,MT) This file is distributed under the terms of the GNU General Public License, see ~ABINIT/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~ABINIT/Infos/contributors.
INPUTS
dmatpuopt= select expression for the density matrix exchmix= mixing factor for local exact-exchange jpawu(ntypat)= value of J llexexch(ntypat)= value of l on which local exact-exchange applies llpawu(ntypat)= value of l on which LDA+U applies ntypat=number of types of atoms in unit cell. pawang <type(pawang_type)>=paw angular mesh and related data %lmax=Maximum value of angular momentum l+1 %gntselect((2*l_max-1)**2,l_max**2,l_max**2)= selection rules for Gaunt coefficients pawprtvol=output printing level for PAW pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data: upawu(ntypat)= value of U use_dmft = 0 no PAW+DMFT, =1 PAW+DMFT useexexch= 0 if no local exact-exchange; 1 if local exact-exchange usepawu= 0 if no LDA+U; 1 if LDA+U
OUTPUT
pawtab <type(pawtab_type)>=paw tabulated data read at start: %ij_proj=nproj*(nproju+1)/2 %klmntomn(4,lmn2_size) = Array giving im, jm ,in, and jn for each klmn=(ilmn,jlmn) %lnproju(nproj)= value of ln for projectors on which paw+u/local exact-exchange acts. %nproju=number of projectors for orbitals on which paw+u/local exact-exchange acts. %phiphjint(pawtabitypat%ij_proj)=Integral of Phi(:,i)*Phi(:,j) for correlated orbitals. %usepawu=0 if no LDA+U; 1 if LDA+U %useexexch=0 if no local exact-exchange; 1 if local exact-exchange === if usepawu>0 %jpawu= value of J %upawu= value of U %vee(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals === if useexexch>0 %fk %vex(2*lpawu+1,:,:,:)=matrix of the screened interaction for correlated orbitals
PARENTS
bethe_salpeter,gstate,m_entropyDMFT,respfn,screening,sigma
CHILDREN
calc_ubare,poisson,simp_gen,wrtout
SOURCE
66 #if defined HAVE_CONFIG_H 67 #include "config.h" 68 #endif 69 70 #include "abi_common.h" 71 72 subroutine pawpuxinit(dmatpuopt,exchmix,f4of2_sla,f6of2_sla,jpawu,llexexch,llpawu,& 73 & ntypat,pawang,pawprtvol,pawrad,pawtab,upawu,use_dmft,useexexch,usepawu,ucrpa) 74 75 use defs_basis 76 use m_errors 77 use m_profiling_abi 78 use m_special_funcs 79 80 use m_pawang, only : pawang_type 81 use m_pawrad, only : pawrad_type, simp_gen, poisson 82 use m_pawtab, only : pawtab_type 83 84 !This section has been created automatically by the script Abilint (TD). 85 !Do not modify the following lines by hand. 86 #undef ABI_FUNC 87 #define ABI_FUNC 'pawpuxinit' 88 use interfaces_14_hidewrite 89 use interfaces_65_paw, except_this_one => pawpuxinit 90 !End of the abilint section 91 92 implicit none 93 94 !Arguments --------------------------------------------- 95 !scalars 96 integer,intent(in) :: dmatpuopt,ntypat,pawprtvol,use_dmft,useexexch,usepawu 97 real(dp),intent(in) :: exchmix 98 type(pawang_type), intent(in) :: pawang 99 integer,optional, intent(in) :: ucrpa 100 !arrays 101 integer,intent(in) :: llexexch(ntypat),llpawu(ntypat) 102 real(dp),intent(in) :: jpawu(ntypat),upawu(ntypat) 103 real(dp),intent(in) :: f4of2_sla(ntypat),f6of2_sla(ntypat) 104 type(pawrad_type),intent(inout) :: pawrad(ntypat) 105 type(pawtab_type),target,intent(inout) :: pawtab(ntypat) 106 107 !Local variables --------------------------------------- 108 !scalars 109 integer :: icount,il,ilmn,isela,iselb,itemp,itypat,iu,j0lmn,jl,jlmn,ju,klm0u 110 integer :: klm0x,klma,klmb,klmn,kln,kln1,kln2,kyc,lcur,lexexch,lkyc,ll,ll1 111 integer :: lmexexch,lmkyc,lmn_size,lmn2_size,lmpawu,lpawu 112 integer :: m1,m11,m2,m21,m3,m31,m4,m41 113 integer :: mesh_size,int_meshsz,mkyc,sz1 114 real(dp) :: ak,f4of2,f6of2,int1,intg 115 character(len=500) :: message 116 !arrays 117 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 118 real(dp),allocatable :: ff(:),fk(:),gg(:) 119 120 ! ************************************************************************* 121 122 DBG_ENTER("COLL") 123 124 if(useexexch==0.and.usepawu==0.and.use_dmft==0) return 125 126 !PAW+U and local exact-exchange restriction 127 if(useexexch>0.and.usepawu>0)then 128 do itypat=1,ntypat 129 if (llpawu(itypat)/=llexexch(itypat).and.llpawu(itypat)/=-1.and.llexexch(itypat)/=-1) then 130 write(message, '(5a,i2,3a)' )& 131 & ' When PAW+U (usepawu>0) and local exact-exchange (exexch>0)',ch10,& 132 & ' are selected together, they must apply on the same',ch10,& 133 & ' angular momentum (lpawu/=lexexch forbidden, here for typat=',itypat,') !',ch10,& 134 & ' Action: correct your input file.' 135 MSG_ERROR(message) 136 end if 137 end do 138 end if 139 140 !Print title 141 if((usepawu>=1.and.usepawu<=4).or.useexexch>0) write(message, '(3a)' ) ch10,ch10," ******************************************" 142 if(usepawu==1) then 143 write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: FLL" 144 else if(usepawu==2) then 145 write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: AMF" 146 else if(usepawu==3) then 147 write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: AMF (alternative)" 148 else if(usepawu==4) then 149 write(message, '(3a)' ) trim(message),ch10," LDA+U Method used: FLL with no spin polarization in the xc functional" 150 end if 151 if(useexexch>0) write(message, '(3a)' ) trim(message),ch10," PAW Local Exact exchange: PBE0" 152 if((usepawu>=1.and.usepawu<=4).or.useexexch>0) & 153 write(message, '(3a)' ) trim(message),ch10," ******************************************" 154 if(use_dmft==0) then 155 call wrtout(ab_out,message,'COLL') 156 call wrtout(std_out, message,'COLL') 157 end if 158 !if(use_dmft>0) then 159 !write(message, '(3a)' ) ch10, " (see DMFT data in log file) " 160 !call wrtout(ab_out,message,'COLL') 161 !endif 162 163 !Loop on atom types 164 do itypat=1,ntypat 165 indlmn => pawtab(itypat)%indlmn 166 lmn_size=pawtab(itypat)%lmn_size 167 lmn2_size=pawtab(itypat)%lmn2_size 168 mesh_size=pawtab(itypat)%mesh_size 169 int_meshsz=pawrad(itypat)%int_meshsz 170 lcur=-1 171 172 ! PAW+U data 173 if (usepawu>0.or.use_dmft>0) then 174 lcur=llpawu(itypat) 175 pawtab(itypat)%lpawu=lcur 176 if(lcur/=-1) then 177 pawtab(itypat)%usepawu=usepawu 178 pawtab(itypat)%upawu=upawu(itypat) 179 pawtab(itypat)%jpawu=jpawu(itypat) 180 pawtab(itypat)%f4of2_sla=f4of2_sla(itypat) 181 pawtab(itypat)%f6of2_sla=f6of2_sla(itypat) 182 else 183 pawtab(itypat)%usepawu=0 184 pawtab(itypat)%upawu=zero 185 pawtab(itypat)%jpawu=zero 186 pawtab(itypat)%f4of2_sla=zero 187 pawtab(itypat)%f6of2_sla=zero 188 end if 189 end if 190 191 ! Local exact-echange data 192 if (useexexch>0) then 193 lcur=llexexch(itypat) 194 pawtab(itypat)%lexexch=lcur 195 pawtab(itypat)%exchmix=exchmix 196 if(pawtab(itypat)%lexexch==-1) pawtab(itypat)%useexexch=0 197 if(pawtab(itypat)%lexexch/=-1) pawtab(itypat)%useexexch=useexexch 198 end if 199 200 ! Select only atoms with +U 201 if(lcur/=-1) then 202 203 ! Compute number of projectors for LDA+U/local exact-exchange/LDA+DMFT 204 icount=count(indlmn(1,1:lmn_size)==lcur) 205 pawtab(itypat)%nproju=icount/(2*lcur+1) 206 if(useexexch>0.and.pawtab(itypat)%nproju>2) then 207 write(message, '(a,a,a)' )& 208 & ' Error on the number of projectors ',ch10,& 209 & ' more than 2 projectors is not allowed for local exact-exchange' 210 MSG_ERROR(message) 211 end if 212 if(pawtab(itypat)%nproju*(2*lcur+1)/=icount) then 213 message = 'pawpuxinit: Error on the number of projectors ' 214 MSG_BUG(message) 215 end if 216 write(message, '(a,a,i4,a,a,i4)' ) ch10,& 217 & ' pawpuxinit : for species ',itypat,ch10,& 218 & ' number of projectors is',pawtab(itypat)%nproju 219 call wrtout(std_out,message,'COLL') 220 221 pawtab(itypat)%ij_proj=pawtab(itypat)%nproju*(pawtab(itypat)%nproju+1)/2 222 223 ! ================================================== 224 ! A-define useful indexes 225 ! -------------------------------------------------- 226 if (allocated(pawtab(itypat)%lnproju)) then 227 ABI_DEALLOCATE(pawtab(itypat)%lnproju) 228 end if 229 ABI_ALLOCATE(pawtab(itypat)%lnproju,(pawtab(itypat)%nproju)) 230 icount=0 231 do ilmn=1,lmn_size 232 if(indlmn(1,ilmn)==lcur) then 233 icount=icount+1 234 itemp=(icount-1)/(2*lcur+1) 235 if (itemp*(2*lcur+1)==icount-1) then 236 pawtab(itypat)%lnproju(itemp+1)=indlmn(5,ilmn) 237 end if 238 end if 239 end do 240 241 if (allocated(pawtab(itypat)%klmntomn)) then 242 ABI_DEALLOCATE(pawtab(itypat)%klmntomn) 243 end if 244 ABI_ALLOCATE(pawtab(itypat)%klmntomn,(4,lmn2_size)) 245 do jlmn=1,lmn_size 246 jl= indlmn(1,jlmn) 247 j0lmn=jlmn*(jlmn-1)/2 248 do ilmn=1,jlmn 249 il= indlmn(1,ilmn) 250 klmn=j0lmn+ilmn 251 pawtab(itypat)%klmntomn(1,klmn)=indlmn(2,ilmn)+il+1 252 pawtab(itypat)%klmntomn(2,klmn)=indlmn(2,jlmn)+jl+1 253 pawtab(itypat)%klmntomn(3,klmn)=indlmn(3,ilmn) 254 pawtab(itypat)%klmntomn(4,klmn)=indlmn(3,jlmn) 255 end do 256 end do 257 258 ! ================================================== 259 ! B-PAW+U: overlap between atomic wavefunctions 260 ! -------------------------------------------------- 261 ! if (usepawu>0) then 262 if(dmatpuopt==1) then 263 write(message, '(4a)' ) ch10,& 264 & ' pawpuxinit : dmatpuopt=1 ',ch10,& 265 & ' PAW+U: dens. mat. constructed by projection on atomic wfn inside PAW augm. region(s)' 266 call wrtout(std_out,message,'COLL') 267 write(message, '(8a)' ) ch10,& 268 & ' pawpuxinit: WARNING: Check that the first partial wave for lpawu:', ch10, & 269 & ' - Is an atomic eigenfunction ',ch10, & 270 & ' - Is normalized ',ch10, & 271 & ' In other cases, choose dmatpuopt=2' 272 call wrtout(std_out,message,'COLL') 273 else if(dmatpuopt==2) then 274 write(message, '(6a)' ) ch10,& 275 & ' pawpuxinit : dmatpuopt=2 ',ch10,& 276 & ' PAW+U: dens. mat. constructed by selecting contribution',ch10,& 277 & ' for each angular momentum to the density (inside PAW augm. region(s))' 278 call wrtout(std_out,message,'COLL') 279 else if(dmatpuopt==3) then 280 write(message, '(a,a,a,a,a,a)' ) ch10,& 281 & ' pawpuxinit : dmatpuopt=3 ',ch10,& 282 & ' PAW+U: dens. mat. constructed by projection on atomic wfn inside PAW augm. region(s)',ch10,& 283 & ' and normalized inside PAW augm. region(s)' 284 call wrtout(std_out,message,'COLL') 285 write(message, '(6a)' ) ch10,& 286 & ' pawpuxinit: WARNING: Check that the first partial wave for lpawu:', ch10, & 287 & ' is an atomic eigenfunction',ch10, & 288 & ' In the other case, choose dmatpuopt=2' 289 call wrtout(std_out,message,'COLL') 290 end if 291 292 ABI_ALLOCATE(ff,(mesh_size)) 293 ff(:)=zero 294 295 if (allocated(pawtab(itypat)%ph0phiint)) then 296 ABI_DEALLOCATE(pawtab(itypat)%ph0phiint) 297 end if 298 if (allocated(pawtab(itypat)%zioneff)) then 299 ABI_DEALLOCATE(pawtab(itypat)%zioneff) 300 end if 301 ABI_ALLOCATE(pawtab(itypat)%ph0phiint,(pawtab(itypat)%nproju)) 302 ABI_ALLOCATE(pawtab(itypat)%zioneff,(pawtab(itypat)%nproju)) 303 304 icount=0 305 do iu=1,pawtab(itypat)%nproju 306 ! write(std_out,*)'DJA iu',iu,' mesh_size',pawtab(itypat)%mesh_size 307 ! do ju=2,pawtab(itypat)%mesh_size 308 ! ff(ju)=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawrad(itypat)%rad(ju) 309 ! write(std_out,fmt='(i5,3e15.5)')ju,pawrad(itypat)%rad(ju),ff(ju),& 310 ! & RadFnH(pawrad(itypat)%rad(ju),4,3,15.0_dp) 311 ! end do 312 ! ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))**2 313 ! call simp_gen(int1,ff,pawrad(itypat)) 314 ! write(std_out,*)'DJA iu',iu,'int1 ',int1 315 ! write(std_out,*)'DJA int1,IRadFnH',int1,IRadFnH(0.0_dp,pawrad(itypat)%rmax,4,3,12) 316 ! Calculation of zioneff 317 ju=pawtab(itypat)%mesh_size-1 318 ak=pawtab(itypat)%phi(ju,pawtab(itypat)%lnproju(iu))/pawtab(itypat)%phi(ju+1,pawtab(itypat)%lnproju(iu)) 319 ak=ak*(pawrad(itypat)%rad(ju+1)/pawrad(itypat)%rad(ju))**(pawtab(itypat)%lpawu-1) 320 pawtab(itypat)%zioneff(iu)=log(ak)/(pawrad(itypat)%rad(ju+1)-pawrad(itypat)%rad(ju)) 321 ! Calculation of ph0phiint 322 ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(1))& 323 & *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu)) 324 call simp_gen(int1,ff,pawrad(itypat)) 325 pawtab(itypat)%ph0phiint(iu)=int1 326 end do 327 if(abs(pawprtvol)>=2) then 328 do icount=1,pawtab(itypat)%nproju 329 write(message, '(a,a,i2,f9.5,a)' ) ch10,& 330 & ' pawpuxinit: icount, ph0phiint(icount)=',icount,pawtab(itypat)%ph0phiint(icount) 331 call wrtout(std_out,message,'COLL') 332 write(message, '(a,f15.5)' ) & 333 & ' pawpuxinit: zioneff=',pawtab(itypat)%zioneff(icount) 334 call wrtout(std_out,message,'COLL') 335 end do 336 write(message, '(a)' ) ch10 337 call wrtout(std_out,message,'COLL') 338 end if 339 340 if (allocated(pawtab(itypat)%phiphjint)) then 341 ABI_DEALLOCATE(pawtab(itypat)%phiphjint) 342 end if 343 ABI_ALLOCATE(pawtab(itypat)%phiphjint,(pawtab(itypat)%ij_proj)) 344 345 icount=0 346 do ju=1,pawtab(itypat)%nproju 347 do iu=1,ju 348 icount=icount+1 349 if ((dmatpuopt==1).and.(useexexch==0)) then 350 pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)*& 351 & pawtab(itypat)%ph0phiint(ju) 352 else if((dmatpuopt==2).or.(useexexch>0)) then 353 ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(iu))& 354 & *pawtab(itypat)%phi(1:mesh_size,pawtab(itypat)%lnproju(ju)) 355 call simp_gen(int1,ff,pawrad(itypat)) 356 pawtab(itypat)%phiphjint(icount)=int1 357 else if((dmatpuopt>=3).and.(useexexch==0)) then 358 pawtab(itypat)%phiphjint(icount)=pawtab(itypat)%ph0phiint(iu)* & 359 & pawtab(itypat)%ph0phiint(ju)/pawtab(itypat)%ph0phiint(1)**(dmatpuopt-2) 360 else 361 write(message, '(3a)' )& 362 & ' PAW+U: dmatpuopt has a wrong value !',ch10,& 363 & ' Action : change value in input file' 364 MSG_ERROR(message) 365 end if 366 end do 367 end do 368 if(pawtab(itypat)%ij_proj/=icount) then 369 message = ' Error in the loop for calculating phiphjint ' 370 MSG_ERROR(message) 371 end if 372 ABI_DEALLOCATE(ff) 373 if(abs(pawprtvol)>=2) then 374 do icount=1,pawtab(itypat)%ij_proj 375 write(message, '(a,a,i2,f9.5,a)' ) ch10,& 376 & ' PAW+U: icount, phiphjint(icount)=',icount,pawtab(itypat)%phiphjint(icount) 377 call wrtout(std_out,message,'COLL') 378 end do 379 end if 380 ! end if 381 382 ! ====================================================================== 383 ! C-PAW+U: Matrix elements of coulomb interaction (see PRB vol.52 5467) 384 ! 1. angular part computed from Gaunt coefficients 385 ! -------------------------------------------------------------------- 386 if (usepawu>0) then 387 lpawu=lcur 388 389 ! a. compute F(k) 390 ! --------------------------------------------- 391 ABI_ALLOCATE(fk,(lpawu+1)) 392 fk(1)=pawtab(itypat)%upawu 393 ! cf Slater Physical Review 165, p 665 (1968) 394 ! write(std_out,*) "f4of2_sla",pawtab(itypat)%f4of2_sla 395 if(lpawu==0) then 396 fk(1)=fk(1) 397 else if(lpawu==1) then 398 fk(2)=pawtab(itypat)%jpawu*5._dp 399 else if(lpawu==2) then 400 ! f4of2=0._dp 401 if(pawtab(itypat)%f4of2_sla<-0.1_dp) then 402 f4of2=0.625_dp 403 pawtab(itypat)%f4of2_sla=f4of2 404 else 405 f4of2=pawtab(itypat)%f4of2_sla 406 end if 407 fk(2)=pawtab(itypat)%jpawu*14._dp/(One+f4of2) 408 fk(3)=fk(2)*f4of2 409 ! if(abs(pawprtvol)>=2) then 410 write(message,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,& 411 & "Slater parameters F^0, F^2, F^4 are",fk(1),fk(2),fk(3) 412 call wrtout(std_out,message,'COLL') 413 ! end if 414 else if(lpawu==3) then 415 f4of2=0.6681_dp 416 f6of2=0.4943_dp 417 if(pawtab(itypat)%f4of2_sla<-0.1_dp) then 418 f4of2=0.6681_dp 419 pawtab(itypat)%f4of2_sla=f4of2 420 else 421 f4of2=pawtab(itypat)%f4of2_sla 422 end if 423 if(pawtab(itypat)%f6of2_sla<-0.1_dp) then 424 f6of2=0.4943_dp 425 pawtab(itypat)%f6of2_sla=f6of2 426 else 427 f6of2=pawtab(itypat)%f6of2_sla 428 end if 429 fk(2)=pawtab(itypat)%jpawu*6435._dp/(286._dp+195._dp*f4of2+250._dp*f6of2) 430 fk(3)=fk(2)*f4of2 431 fk(4)=fk(2)*f6of2 432 write(std_out,'(a,3x,a,f9.4,f9.4,f9.4,f9.4)') ch10,& 433 & "Slater parameters F^0, F^2, F^4, F^6 are",fk(1),fk(2),fk(3),fk(4) 434 else 435 write(message, '(a,i0,2a)' )& 436 & ' lpawu=',lpawu,ch10,& 437 & ' lpawu not equal to 0 ,1 ,2 or 3 is not allowed' 438 MSG_ERROR(message) 439 end if 440 441 ! b. Compute ak and vee. 442 ! --------------------------------------------- 443 if (allocated(pawtab(itypat)%vee)) then 444 ABI_DEALLOCATE(pawtab(itypat)%vee) 445 end if 446 sz1=2*lpawu+1 447 ABI_ALLOCATE(pawtab(itypat)%vee,(sz1,sz1,sz1,sz1)) 448 pawtab(itypat)%vee=zero 449 lmpawu=(lpawu-1)**2+2*(lpawu-1)+1 ! number of m value below correlated orbitals 450 klm0u=lmpawu*(lmpawu+1)/2 ! value of klmn just below correlated orbitals 451 ! --------- 4 loops for interaction matrix 452 do m1=-lpawu,lpawu 453 m11=m1+lpawu+1 454 do m2=-lpawu,m1 455 m21=m2+lpawu+1 456 ! klma= number of pair before correlated orbitals + 457 ! number of pair for m1 lower than correlated orbitals 458 ! (m1+lpawu+1)*(lpawu-1) + number of pairs for correlated orbitals 459 ! before (m1,m2) + number of pair for m2 lower than current value 460 klma=klm0u+m11*lmpawu+(m11-1)*m11/2+m21 461 do m3=-lpawu,lpawu 462 m31=m3+lpawu+1 463 do m4=-lpawu,m3 464 m41=m4+lpawu+1 465 klmb=klm0u+m31*lmpawu+(m31-1)*m31/2+m41 466 ! --------- loop on k=1,2,3 (4 if f orbitals) 467 do kyc=1,2*lpawu+1,2 468 lkyc=kyc-1 469 lmkyc=(lkyc+1)*(lkyc)+1 470 ak=zero 471 do mkyc=-lkyc,lkyc,1 472 isela=pawang%gntselect(lmkyc+mkyc,klma) 473 iselb=pawang%gntselect(lmkyc+mkyc,klmb) 474 if (isela>0.and.iselb>0) ak=ak +pawang%realgnt(isela)*pawang%realgnt(iselb) 475 end do 476 ! ----- end loop on k=1,2,3 (4 if f orbitals) 477 ak=ak/(two*dble(lkyc)+one) 478 pawtab(itypat)%vee(m11,m31,m21,m41)=ak*fk(lkyc/2+1)+pawtab(itypat)%vee(m11,m31,m21,m41) 479 end do !kyc 480 pawtab(itypat)%vee(m11,m31,m21,m41)=pawtab(itypat)%vee(m11,m31,m21,m41)*four_pi 481 pawtab(itypat)%vee(m21,m31,m11,m41)=pawtab(itypat)%vee(m11,m31,m21,m41) 482 pawtab(itypat)%vee(m11,m41,m21,m31)=pawtab(itypat)%vee(m11,m31,m21,m41) 483 pawtab(itypat)%vee(m21,m41,m11,m31)=pawtab(itypat)%vee(m11,m31,m21,m41) 484 485 !!!! pawtab(itypat)%vee(m11,m31,m21,m41)= <m11 m31| vee| m21 m41 > 486 end do 487 end do 488 end do 489 end do 490 ABI_DEALLOCATE(fk) 491 ! testu=0 492 ! write(std_out,*) " Matrix of interaction vee(m1,m2,m1,m2)" 493 ! do m1=1,2*lpawu+1 494 ! write(std_out,'(2x,14(f12.6,2x))') (pawtab(itypat)%vee(m1,m2,m1,m2),m2=1,2*lpawu+1) 495 ! do m2=1,2*lpawu+1 496 ! testu=testu+ pawtab(itypat)%vee(m1,m2,m1,m2) 497 ! enddo 498 ! enddo 499 ! testu=testu/((two*lpawu+one)**2) 500 ! write(std_out,*) "------------------------" 501 ! write(std_out,'(a,f12.6)') " U=", testu 502 ! write(std_out,*) "------------------------" 503 ! write(std_out,*) " Matrix of interaction vee(m1,m2,m1,m2)-vee(m1,m2,m2,m1)" 504 ! do m1=1,2*lpawu+1 505 ! write(std_out,'(2x,14(f12.6,2x))') ((pawtab(itypat)%vee(m1,m2,m1,m2)-pawtab(itypat)%vee(m1,m2,m2,m1)),m2=1,2*lpawu+1) 506 ! do m2=1,2*lpawu+1 507 ! if(m1/=m2) testumj=testumj+ pawtab(itypat)%vee(m1,m2,m1,m2)-pawtab(itypat)%vee(m1,m2,m2,m1) 508 ! enddo 509 ! enddo 510 ! testumj=testumj/((two*lpawu)*(two*lpawu+one)) 511 ! write(std_out,*) "------------------------" 512 ! write(std_out,'(a,f12.6)') " U-J=", testumj 513 ! write(std_out,*) "------------------------" 514 ! write(std_out,*) "------------------------" 515 ! write(std_out,'(a,f12.6)') " J=", testu-testumj 516 ! write(std_out,*) "------------------------" 517 end if ! usepawu 518 519 ! ====================================================================== 520 ! D-Local ex-exchange: Matrix elements of coulomb interaction and Fk 521 ! ---------------------------------------------------------------------- 522 if (useexexch>0) then 523 lexexch=lcur 524 525 ! a. compute F(k) 526 ! --------------------------------------------- 527 if (allocated(pawtab(itypat)%fk)) then 528 ABI_DEALLOCATE(pawtab(itypat)%fk) 529 end if 530 ABI_ALLOCATE(pawtab(itypat)%fk,(6,4)) 531 pawtab(itypat)%fk=zero 532 ABI_ALLOCATE(ff,(mesh_size)) 533 ABI_ALLOCATE(gg,(mesh_size)) 534 ff(:)=zero;gg(:)=zero 535 kln=(pawtab(itypat)%lnproju(1)*( pawtab(itypat)%lnproju(1)+1)/2) 536 do ll=1,lexexch+1 537 ll1=2*ll-2 538 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 539 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln) 540 call poisson(ff,ll1,pawrad(itypat),gg) 541 ff(1)=zero 542 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln)*gg(2:mesh_size))& 543 & /pawrad(itypat)%rad(2:mesh_size) 544 call simp_gen(intg,ff,pawrad(itypat)) 545 pawtab(itypat)%fk(1,ll)=intg*(two*ll1+one) 546 end do 547 if (pawtab(itypat)%nproju==2) then 548 kln1=kln+pawtab(itypat)%lnproju(1) 549 kln2=kln1+1 550 do ll=1,lexexch+1 551 ll1=2*ll-2 552 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 553 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln1) 554 call poisson(ff,ll1,pawrad(itypat),gg) 555 ff(1)=zero 556 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln1)*gg(2:mesh_size))& 557 & /pawrad(itypat)%rad(2:mesh_size) 558 call simp_gen(intg,ff,pawrad(itypat)) 559 pawtab(itypat)%fk(2,ll)=intg*(two*ll1+one) 560 end do 561 do ll=1,lexexch+1 562 ll1=2*ll-2 563 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 564 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln2) 565 call poisson(ff,ll1,pawrad(itypat),gg) 566 ff(1)=zero 567 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))& 568 & /pawrad(itypat)%rad(2:mesh_size) 569 call simp_gen(intg,ff,pawrad(itypat)) 570 pawtab(itypat)%fk(3,ll)=intg*(two*ll1+one) 571 end do 572 do ll=1,lexexch+1 573 ll1=2*ll-2 574 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 575 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln) 576 call poisson(ff,ll1,pawrad(itypat),gg) 577 ff(1)=zero 578 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln1)*gg(2:mesh_size))& 579 & /pawrad(itypat)%rad(2:mesh_size) 580 call simp_gen(intg,ff,pawrad(itypat)) 581 pawtab(itypat)%fk(4,ll)=intg*(two*ll1+one) 582 end do 583 do ll=1,lexexch+1 584 ll1=2*ll-2 585 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 586 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln) 587 call poisson(ff,ll1,pawrad(itypat),gg) 588 ff(1)=zero 589 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))& 590 & /pawrad(itypat)%rad(2:mesh_size) 591 call simp_gen(intg,ff,pawrad(itypat)) 592 pawtab(itypat)%fk(5,ll)=intg*(two*ll1+one) 593 end do 594 do ll=1,lexexch+1 595 ll1=2*ll-2 596 if (int_meshsz<mesh_size) ff(int_meshsz+1:mesh_size)=zero 597 ff(1:int_meshsz)=pawtab(itypat)%phiphj(1:int_meshsz,kln1) 598 call poisson(ff,ll1,pawrad(itypat),gg) 599 ff(1)=zero 600 ff(2:mesh_size)=(pawtab(itypat)%phiphj(2:mesh_size,kln2)*gg(2:mesh_size))& 601 & /pawrad(itypat)%rad(2:mesh_size) 602 call simp_gen(intg,ff,pawrad(itypat)) 603 pawtab(itypat)%fk(6,ll)=intg*(two*ll1+one) 604 end do 605 f4of2=0.6681_dp 606 f6of2=0.4943_dp 607 end if 608 ABI_DEALLOCATE(ff) 609 ABI_DEALLOCATE(gg) 610 611 ! b. Compute vex. 612 ! --------------------------------------------- 613 if (allocated(pawtab(itypat)%vex)) then 614 ABI_DEALLOCATE(pawtab(itypat)%vex) 615 end if 616 sz1=2*lexexch+1 617 ABI_ALLOCATE(pawtab(itypat)%vex,(sz1,sz1,sz1,sz1,4)) 618 pawtab(itypat)%vex=zero 619 lmexexch=(lexexch-1)**2+2*(lexexch-1)+1 ! number of m value below correlated orbitals 620 klm0x=lmexexch*(lmexexch+1)/2 ! value of klmn just below correlated orbitals 621 ! --------- 4 loops for interaction matrix 622 do m1=-lexexch,lexexch 623 m11=m1+lexexch+1 624 do m2=-lexexch,m1 625 m21=m2+lexexch+1 626 ! klma= number of pair before correlated orbitals + 627 ! number of pair for m1 lower than correlated orbitals 628 ! (m1+lexexch+1)*(lexexch-1) + number of pairs for correlated orbitals 629 ! before (m1,m2) + number of pair for m2 lower than current value 630 klma=klm0x+m11*lmexexch+(m11-1)*m11/2+m21 631 do m3=-lexexch,lexexch 632 m31=m3+lexexch+1 633 do m4=-lexexch,m3 634 m41=m4+lexexch+1 635 klmb=klm0x+m31*lmexexch+(m31-1)*m31/2+m41 636 ! --------- loop on k=1,2,3 (4 if f orbitals) 637 do kyc=1,2*lexexch+1,2 638 lkyc=kyc-1 639 ll=(kyc+1)/2 640 lmkyc=(lkyc+1)*(lkyc)+1 641 ak=zero 642 do mkyc=-lkyc,lkyc,1 643 isela=pawang%gntselect(lmkyc+mkyc,klma) 644 iselb=pawang%gntselect(lmkyc+mkyc,klmb) 645 if (isela>0.and.iselb>0) ak=ak +pawang%realgnt(isela)*pawang%realgnt(iselb) 646 end do 647 ! ----- end loop on k=1,2,3 (4 if f orbitals) 648 pawtab(itypat)%vex(m11,m31,m21,m41,ll)=ak/(two*dble(lkyc)+one) 649 end do !kyc 650 do ll=1,4 651 pawtab(itypat)%vex(m11,m31,m21,m41,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll)*four_pi 652 pawtab(itypat)%vex(m21,m31,m11,m41,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll) 653 pawtab(itypat)%vex(m11,m41,m21,m31,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll) 654 pawtab(itypat)%vex(m21,m41,m11,m31,ll)=pawtab(itypat)%vex(m11,m31,m21,m41,ll) 655 end do 656 end do 657 end do 658 end do 659 end do 660 661 end if !useexexch>0 662 663 if (present(ucrpa)) then 664 if (ucrpa>=1) then 665 call calc_ubare(itypat,lcur,pawang,pawrad(itypat),pawtab(itypat)) 666 call calc_ubare(itypat,lcur,pawang,pawrad(itypat),pawtab(itypat),pawtab(itypat)%rpaw) 667 end if 668 end if 669 end if !lcur/=-1 670 end do !end loop on typat 671 672 DBG_EXIT("COLL") 673 674 end subroutine pawpuxinit