TABLE OF CONTENTS
- ABINIT/datafordmft
- ABINIT/m_datafordmft
- m_datafordmft/compute_levels
- m_datafordmft/hybridization_asymptotic_coefficient
- m_datafordmft/psichi_check
- m_datafordmft/psichi_print
- m_datafordmft/psichi_renormalization
- psichi_renormalization/normalizepsichi
ABINIT/datafordmft [ Functions ]
NAME
datafordmft
FUNCTION
Compute psichi (and print some data for check)
INPUTS
cryst_struc <type(crystal_t)>=crystal structure data -gprimd(3,3)=dimensional reciprocal space primitive translations -indsym(4,nsym,natom)=indirect indexing array for atom labels -symrec(3,3,nsym)=symmetry operations in reciprocal space - nsym= number of symetry operations cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector dimcprj(natom) = dimension for cprj dtset <type(dataset_type)>=all input variables for this dataset eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) fermie= Fermi energy dft_occup <type(oper_type)> = occupations in the correlated orbitals in DFT mband=maximum number of bands mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization nkpt=number of k points. my_nspinor=number of spinorial components of the wavefunctions (on current proc) nsppol=1 for unpolarized, 2 for spin-polarized occ(mband*nkpt*nsppol) = occupancies of KS states. paw_dmft <type(paw_dmft_type)>= paw+dmft related data paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang)>=paw angular mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials unpaw = file number for cprj
OUTPUT
paw_dmft%psichi(nsppol,nkpt,mband,my_nspinor,dtset%natom,(2*maxlpawu+1))): projections <Psi|chi> paw_dmft%eigen(paw_dmft%nsppol,paw_dmft%nkpt,paw_dmft%mband)
SIDE EFFECTS
(only writing, printing)
NOTES
SOURCE
107 subroutine datafordmft(cryst_struc,cprj,dimcprj,dtset,eigen,fermie,& 108 & dft_occup,mband,mband_cprj,mkmem,mpi_enreg,nkpt,my_nspinor,nsppol,occ,& 109 & paw_dmft,paw_ij,pawang,pawtab,psps,usecprj,unpaw,nbandkss) 110 111 !Arguments ------------------------------------ 112 !scalars 113 integer,intent(in) :: mband,mband_cprj,mkmem 114 integer,intent(in) :: nkpt,my_nspinor,nsppol 115 integer,intent(in) :: unpaw,usecprj 116 integer, optional, intent(in) :: nbandkss 117 real(dp),intent(in) :: fermie 118 type(MPI_type),intent(in) :: mpi_enreg 119 type(dataset_type),intent(in) :: dtset 120 type(oper_type), intent(inout) :: dft_occup !vz_i 121 type(pawang_type),intent(in) :: pawang 122 type(pseudopotential_type),intent(in) :: psps 123 type(crystal_t),intent(in) :: cryst_struc 124 !arrays 125 integer, intent(in) :: dimcprj(cryst_struc%natom) 126 real(dp),intent(in) :: eigen(mband*nkpt*nsppol) 127 real(dp),intent(in) :: occ(mband*nkpt*nsppol) 128 type(paw_ij_type),intent(in) :: paw_ij(:) 129 ! type(pawcprj_type) :: cprj(cryst_struc%natom,my_nspinor*mband*mkmem*nsppol) 130 type(pawcprj_type), intent(in) :: cprj(cryst_struc%natom,my_nspinor*mband_cprj*mkmem*nsppol*usecprj) 131 type(paw_dmft_type), intent(inout) :: paw_dmft 132 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 133 !Local variables------------------------------- 134 !scalars 135 integer :: band_index,dimpsichi,facpara 136 integer :: iat,iatom,ib,iband,ibandc,ibg,icat,icount_proj_ilmn,idijeff,ierr,ierrr,ikpt 137 integer :: ilmn,im,im1,iorder_cprj,ispinor,ispinor1,isppol,itypat,ilmn1 138 integer :: jj1,ldim,lmn_size,lpawu 139 integer :: m1,m1_x2my2d,m1_t2g,m1_t2g_mod,maxnproju,me,natom,nband_k,nband_k_cprj 140 integer :: nbandi,nbandf,nnn,nprocband,nsploop,option,opt_renorm,spaceComm,unt 141 real(dp) :: ph0phiint_used 142 character(len=500) :: message 143 !arrays 144 real(dp) :: chinorm 145 complex(dpc), allocatable :: buffer1(:) 146 type(matlu_type), allocatable :: loc_occ_check(:) 147 type(matlu_type), allocatable :: loc_norm_check(:) 148 type(matlu_type), allocatable :: xocc_check(:) 149 type(matlu_type), allocatable :: xnorm_check(:) 150 type(matlu_type), allocatable :: matlu_temp(:) 151 logical :: lprojchi,t2g,x2my2d 152 integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/)) 153 type(pawcprj_type),allocatable :: cwaveprj(:,:) 154 !************************************************************************ 155 156 !DBG_ENTER("COLL") 157 !Fake test to keep fermie as argument. REMOVE IT AS SOON AS POSSIBLE ... 158 if(fermie>huge(zero))chinorm=zero 159 160 facpara=1 !mpi_enreg%nproc 161 if(abs(dtset%pawprtvol)>=3) then 162 write(message,*) " number of k-points used is nkpt=nkpt ", nkpt 163 call wrtout(std_out, message,'COLL') 164 write(message,*) " warning: parallelised version ", nkpt 165 call wrtout(std_out, message,'COLL') 166 write(message,*) " weights k-points used is wtk=wtk" 167 call wrtout(std_out, message,'COLL') 168 end if 169 170 if(usecprj==0) then 171 write(message,*) " usecprj=0 : BUG in datafordmft",usecprj 172 ABI_BUG(message) 173 end if 174 175 if(my_nspinor/=dtset%nspinor) then 176 write(message,*) " my_nspinor=/dtset%nspinor, datafordmft not working in this case",& 177 & my_nspinor,dtset%nspinor 178 ABI_ERROR(message) 179 end if 180 181 182 !do ib=1,my_nspinor*mband_cprj*mkmem*nsppol*usecprj 183 !write(std_out,'(a,i6,3e16.7)') "cprj",ib,cprj(1,ib)%cp(1,19),cprj(1,ib)%cp(2,19),cprj(1,ib)%cp(1,19)**2+cprj(1,ib)%cp(2,19)**2 184 !enddo 185 186 !----------------------------------- MPI------------------------------------- 187 188 !Init parallelism 189 spaceComm=mpi_enreg%comm_cell 190 if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt 191 me=mpi_enreg%me_kpt 192 iorder_cprj=0 193 ABI_CHECK(dtset%mkmem/=0,"mkmem==0 not supported anymore!") 194 !todo_ab: extract cprj from file unpaw in the following.. 195 !call abi_abort('COLL') 196 197 !----------------------------------- MPI------------------------------------- 198 199 nbandi=paw_dmft%dmftbandi 200 nbandf=paw_dmft%dmftbandf 201 lprojchi=.false. 202 lprojchi=.true. 203 t2g=(paw_dmft%dmftqmc_t2g==1) 204 x2my2d=(paw_dmft%dmftqmc_x2my2d==1) 205 natom=cryst_struc%natom 206 207 !if(mpi_enreg%me==0) write(7886,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc 208 !if(mpi_enreg%me==1) write(7887,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc 209 !if(mpi_enreg%me==2) write(7888,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc 210 write(message,'(2a)') ch10,& 211 & ' == Prepare data for DMFT calculation ' 212 call wrtout(std_out,message,'COLL') 213 if(abs(dtset%pawprtvol)>=3) then 214 write(message, '(a,a)' ) ch10,& 215 & '---------------------------------------------------------------' 216 ! call wrtout(ab_out,message,'COLL') 217 ! call wrtout(std_out, message,'COLL') 218 call wrtout(std_out, message,'COLL') 219 write(message, '(a,a,a,a,a,a,a,a,a,a,a,a)' ) ch10,& 220 & ' Print useful data (as a check)',ch10,& 221 & ' - Overlap of KS wfc with atomic orbital inside sphere',ch10,& 222 & ' - Eigenvalues',ch10,& 223 & ' - Weights of k-points',ch10,& 224 & ' - Number of spins ',ch10,& 225 & ' - Number of states' 226 ! call wrtout(ab_out,message,'COLL') 227 ! call wrtout(std_out, message,'COLL') 228 call wrtout(std_out, message,'COLL') 229 write(message, '(a,a)' ) ch10,& 230 & '---------------------------------------------------------------' 231 call wrtout(std_out, message, 'COLL') 232 end if 233 if(dtset%nstep==0.and.dtset%nbandkss==0) then 234 message = 'nstep should be greater than 1' 235 ABI_BUG(message) 236 end if 237 238 !********************* Max Values for U terms. 239 !maxlpawu=0 240 maxnproju=0 241 do iatom=1,natom 242 if(pawtab(dtset%typat(iatom))%lpawu.ne.-1 .and. pawtab(dtset%typat(iatom))%nproju.gt.maxnproju)& 243 & maxnproju=pawtab(dtset%typat(iatom))%nproju 244 end do 245 !***************** in forlb.eig 246 if(me.eq.0.and.abs(dtset%pawprtvol)>=3) then 247 if (open_file('forlb.eig',message,newunit=unt,form='formatted',status='unknown') /= 0) then 248 ABI_ERROR(message) 249 end if 250 rewind(unt) 251 write(unt,*) "Number of bands, spins, and k-point; and spin-orbit flag" 252 write(unt,*) mband,nsppol,nkpt,my_nspinor,nbandi,nbandf 253 write(unt,*) " For each k-point, eigenvalues for each band" 254 write(unt,*) (dtset%wtk(ikpt)*facpara,ikpt=1,nkpt) 255 band_index=0 256 do isppol=1,nsppol 257 write(unt,*) " For spin" 258 write(unt,*) isppol 259 do ikpt=1,nkpt 260 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 261 ibandc=0 262 write(unt,*) " For k-point" 263 write(unt,*) ikpt 264 do iband=1,mband 265 if(paw_dmft%band_in(iband)) then 266 ibandc=ibandc+1 267 write(unt, '(2i6,4x,f20.15)' ) ibandc,ikpt,eigen(iband+band_index)*2.d0 268 end if 269 end do 270 band_index=band_index+nband_k 271 end do 272 end do 273 close(unt) 274 end if ! proc=me 275 276 !== put eigen into eigen_dft 277 band_index=0 278 do isppol=1,nsppol 279 do ikpt=1,nkpt 280 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 281 ibandc=0 282 do iband=1,mband 283 if(paw_dmft%band_in(iband)) then 284 ibandc=ibandc+1 285 paw_dmft%eigen_dft(isppol,ikpt,ibandc)=eigen(iband+band_index) ! in Ha 286 ! paw_dmft%eigen_dft(isppol,ikpt,ibandc)=fermie 287 end if 288 end do 289 band_index=band_index+nband_k 290 end do 291 end do 292 293 if(abs(dtset%pawprtvol)>=3) then 294 write(message, '(a,a)' ) ch10,& 295 & ' datafordmft : eigenvalues written' 296 call wrtout(std_out, message,'COLL') 297 end if 298 !========================================================================== 299 !***************** Compute <Psi|Chi>=\sum_{proja} <Psi|P_a><phi_a|Chi> 300 !========================================================================== 301 !write(std_out,*) "size(cprj,dim=1)",size(cprj,dim=1),size(cprj,dim=2),dtset%mband,dtset%mkmem,dtset%nkpt 302 303 !Allocate temporary cwaveprj storage 304 ABI_MALLOC(cwaveprj,(natom,my_nspinor)) 305 !write(std_out,*) "before alloc cprj" 306 !write(std_out,*) size(cwaveprj,dim=1),size(cwaveprj,dim=2),size(dimcprj,dim=1) 307 308 call pawcprj_alloc(cwaveprj,0,dimcprj) 309 !write(std_out,*) "after alloc cprj" 310 311 nprocband=(mband/mband_cprj) 312 ibg=0 313 do isppol=1,nsppol 314 do ikpt=1,nkpt 315 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 316 nband_k_cprj=nband_k/nprocband 317 318 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 319 320 ! write(2011,*) ikpt 321 ! ibsp=ibg 322 ibandc=0 323 ! LOOP OVER BANDS 324 ib=0 325 do iband=1,nband_k 326 327 if(paw_dmft%band_in(iband)) ibandc=ibandc+1 328 329 ! Parallelization: treat only some bands 330 if (dtset%paral_kgb==1) then 331 if (mod((iband-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle 332 else 333 if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle 334 end if 335 ib=ib+1 336 337 do ispinor=1,my_nspinor 338 ! ibsp =~ (ikpt-1)*nband*my_nspinor+iband 339 ! ibsp=ibsp+1 340 icat=1 341 ! write(std_out,*) isppol,ikpt,iband,ispinor 342 iat=0 ! to declare 343 do itypat=1,dtset%ntypat 344 lmn_size=pawtab(itypat)%lmn_size 345 ! write(std_out,*) isppol,ikpt,iband,ispinor 346 do iatom=icat,icat+cryst_struc%nattyp(itypat)-1 347 lpawu=pawtab(itypat)%lpawu 348 ! ---------- t2g case 349 if(paw_dmft%dmftqmc_t2g==1.and.lpawu/=-1) then 350 if(lpawu==2) then 351 ! lpawu==2 must be chosen in input and thus in 352 ! pawtab. On the contrary, paw_dmft now has 353 ! lpawu=1 354 m1_t2g=0 ! index for psichi which has a dimension 3 355 else 356 write(message,'(a,a,i4,i4,2a)') ch10,& 357 & ' Wrong use of dmftqmc_t2g',paw_dmft%dmftqmc_t2g,lpawu,ch10,& 358 & ' Action: desactivate qmftqmc_t2g or use lpawu=2' 359 ABI_ERROR(message) 360 end if 361 end if 362 ! ---------- t2g case 363 ! ---------- x2my2d case 364 if(paw_dmft%dmftqmc_x2my2d==1.and.lpawu/=-1) then 365 if(lpawu==2) then 366 ! lpawu==2 must be chosen in input and thus in 367 ! pawtab. On the contrary, paw_dmft now has 368 ! lpawu=1 369 m1_x2my2d=0 ! index for psichi which has a dimension 1 370 else 371 write(message,'(a,a,i4,i4,2a)') ch10,& 372 & ' Wrong use of dmftqmc_x2my2d',paw_dmft%dmftqmc_x2my2d,lpawu,ch10,& 373 & ' Action: desactivate dmftqmc_x2my2d or use lpawu=2' 374 ABI_ERROR(message) 375 end if 376 end if 377 ! ---------- x2my2d case 378 ! if(isppol==2) write(std_out,*) "ee",size(cprj(iatom,ibsp)%cp(:,:)) 379 380 iat=iat+1 381 jj1=0 382 if(lpawu.ne.-1) then 383 384 call pawcprj_get(cryst_struc%atindx1,cwaveprj,cprj,natom,ib,ibg,ikpt,& 385 & iorder_cprj,isppol,mband_cprj,mkmem,natom,1,nband_k_cprj,& 386 & my_nspinor,nsppol,unpaw,mpicomm=mpi_enreg%comm_kpt,& 387 & proc_distrb=mpi_enreg%proc_distrb) 388 ! write(std_out,*) "M-2",isppol,ikpt,iband,iatom,& 389 ! & (cwaveprj(iatom,ispinor)%cp(1,13)**2+cwaveprj(iatom,ispinor)%cp(2,13)**2) ,ibsp 390 391 ! chinorm=(pawtab(itypat)%phiphjint(1)) 392 chinorm=1.d0 393 ! write(std_out,*) isppol,ikpt,iband,ispinor,iat 394 do ilmn=1,lmn_size 395 ! write(std_out,*) ilmn 396 ! ------------ Select l=lpawu. 397 if (psps%indlmn(1,ilmn,itypat)==lpawu) then 398 ! ------------ Check that the band is choosen (within nbandi and nbandf) 399 if(paw_dmft%band_in(iband)) then 400 ! if(ilmn==13) then 401 ! write(std_out,*) "M-2c",isppol,ikpt,iband,iatom 402 ! write(std_out,*) "M-2b",isppol,ikpt,ibandc,& 403 ! & (cwaveprj(iatom,ispinor)%cp(1,13)**2+cwaveprj(iatom,ispinor)%cp(2,13)**2) 404 ! endif 405 ! write(std_out,*) "inside paw_dmft%band_in",iband 406 jj1=jj1+1 407 if(jj1>pawtab(itypat)%nproju*(2*lpawu+1)) then 408 write(message,'(a,a,a,a)') ch10,& 409 & ' jj1 is not correct in datafordmft',ch10,& 410 & ' Action: CONTACT Abinit group' 411 ABI_BUG(message) 412 end if ! jj1 413 icount_proj_ilmn=psps%indlmn(3,ilmn,itypat) ! n: nb of projector 414 ! write(std_out,*) "icount_proj_ilmn",icount_proj_ilmn 415 m1=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1 416 ! ---- if lprochi=true do the sum over every projectors 417 ! ---- if lprochi=false: . do the sum over only ground state projector 418 ! ---- . and take ph0phiint(icount_proj_ilmn)=1 419 420 ! call pawcprj_get(cryst_struc%atindx1,cwaveprj,cprj,natom,& 421 ! & ib,ibg,ikpt,iorder_cprj,isppol,mband_cprj,mkmem,& 422 ! & natom,1,nband_k_cprj,my_nspinor,nsppol,unpaw,& 423 ! & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 424 425 if(lprojchi.or.icount_proj_ilmn==1) then 426 if(.not.lprojchi) ph0phiint_used=one 427 if(lprojchi) ph0phiint_used=pawtab(itypat)%ph0phiint(icount_proj_ilmn) 428 if(paw_dmft%dmftqmc_t2g==1) then ! t2g case 429 430 if(m1==1.or.m1==2.or.m1==4) then ! t2g case 431 m1_t2g=m1_t2g+1 ! t2g case1 432 m1_t2g_mod=mod(m1_t2g-1,3)+1 433 ! write(std_out,*) "M0",isppol,ikpt,iband,ilmn,cprj(iatom,ibsp)%cp(1,ilmn) 434 ! write(std_out,*) "M1",m1,m1_t2g,iband,ilmn,icount_proj_ilmn 435 paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g_mod)=& ! t2g case 436 & paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g_mod)+& ! t2g case 437 ! & cmplx(cprj(iatom,ibsp)%cp(1,ilmn)*ph0phiint_used,& ! t2g case 438 ! & cprj(iatom,ibsp)%cp(2,ilmn)*ph0phiint_used,kind=dp) ! t2g case 439 & cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn)*ph0phiint_used,& ! t2g case 440 & cwaveprj(iatom,ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp) ! t2g case 441 end if !t2g case 442 else if(paw_dmft%dmftqmc_x2my2d==1) then ! x2my2d case 443 if(m1==5) then ! x2my2d case 444 m1_x2my2d=1 ! x2my2d case1 445 paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_x2my2d)=& ! x2my2d case 446 & paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_x2my2d)+& ! x2my2d case 447 & cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn)*ph0phiint_used,& ! x2my2d case 448 & cwaveprj(iatom,ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp) ! x2my2d case 449 end if !x2my2d case 450 else 451 paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)=& 452 & paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)+& 453 ! & cmplx(cprj(iatom,ibsp)%cp(1,ilmn)*ph0phiint_used,& 454 ! & cprj(iatom,ibsp)%cp(2,ilmn)*ph0phiint_used,kind=dp) 455 & cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn)*ph0phiint_used,& 456 & cwaveprj(iatom,ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp) 457 458 ! if(ibandc==3.and.iat==1.and.m1==1) then 459 ! write(std_out,'(a,3i5)') "psichi integers",isppol,ikpt,ispinor 460 ! write(std_out,'(a,2i5,2e16.7)') "psichi IB3 iAT1 IM1",ilmn,icount_proj_ilmn,& 461 ! & real(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)), imag(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)) 462 ! write(std_out,'(a,2i5,2e16.7)') "cwaveprj IB3 iAT1 IM1",ilmn,icount_proj_ilmn,cwaveprj(iatom,ispinor)%cp(1,ilmn) & 463 ! & , cwaveprj(iatom,ispinor)%cp(2,ilmn) 464 ! endif 465 end if 466 end if ! lprojchi=.true. (always) 467 end if ! nband belongs to paw_dmft%band_in 468 end if ! L=lpawu 469 end do !ilmn : loop over L,M,N 470 end if ! If lpawu/=-1 471 end do ! iatom 472 icat=icat+cryst_struc%nattyp(itypat) 473 end do ! itypat 474 end do ! ispinor 475 end do !iband 476 ibg=ibg+nband_k_cprj*my_nspinor 477 ! bdtot_index=bdtot_index+nband_k ! useless (only for occ) 478 end do !ikpt 479 end do ! isppol 480 !do isppol=1,nsppol 481 !do ikpt=1,nkpt 482 !do ispinor=1,my_nspinor 483 !write(std_out,*) "psichi integers",isppol,ikpt,ispinor 484 !write(std_out,*) "psichi IB3 iAT1 IM1",& 485 !& real(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)), imag(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)) 486 ! 487 !enddo 488 !enddo 489 !enddo 490 !call abi_abort('COLL') 491 if(abs(dtset%pawprtvol)>=3) then 492 write(message,*) "chinorm used here =",chinorm 493 call wrtout(std_out, message,'COLL') 494 end if 495 496 !deallocate temporary cwaveprj/cprj storage 497 call pawcprj_free(cwaveprj) 498 ABI_FREE(cwaveprj) 499 500 !========================================================================== 501 !********************* Gather information for MPI before printing 502 !========================================================================== 503 504 dimpsichi=2*nsppol*nkpt*mband*my_nspinor*natom*(2*paw_dmft%maxlpawu+1) 505 ABI_MALLOC(buffer1,(dimpsichi)) 506 buffer1 = zero 507 nnn=0 508 !write(176,*) "beg",psichi 509 do isppol=1,nsppol 510 do ikpt=1,nkpt 511 do ibandc=1,paw_dmft%mbandc 512 do ispinor=1,my_nspinor 513 do iat=1,natom 514 do m1=1,2*paw_dmft%maxlpawu+1 515 ! do m=1,2 516 nnn=nnn+1 517 buffer1(nnn)=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1) 518 ! enddo 519 end do 520 end do 521 end do 522 end do 523 end do 524 end do 525 call xmpi_barrier(spaceComm) 526 call xmpi_sum(buffer1,spaceComm ,ierr) 527 if (dtset%paral_kgb==1.and.nprocband>1) then 528 call xmpi_sum(buffer1,mpi_enreg%comm_band,ierr) !Build sum over band processors 529 end if 530 call xmpi_barrier(spaceComm) 531 nnn=0 532 do isppol=1,nsppol 533 do ikpt=1,nkpt 534 do ibandc=1,paw_dmft%mbandc 535 do ispinor=1,my_nspinor 536 do iat=1,natom 537 do m1=1,2*paw_dmft%maxlpawu+1 538 ! do m=1,2 539 nnn=nnn+1 540 paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)=buffer1(nnn) 541 ! enddo 542 ! if(ibandc==1) then 543 ! write(6,*) "m1",m1 544 ! write(6,*) "psichi datafordmft",paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1) 545 ! endif 546 end do 547 end do 548 end do 549 end do 550 end do 551 end do 552 ABI_FREE(buffer1) 553 !write(177,*) "end",psichi 554 555 !do isppol=1,nsppol 556 !do ikpt=1,nkpt 557 !do ibandc=1,paw_dmft%mbandc 558 !do ispinor=1,my_nspinor 559 !write(std_out,*) "psichigather",isppol,ikpt,ibandc,& 560 !& real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,1,1))**2+& 561 !& imag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,1,1))**2 562 ! 563 !enddo 564 !enddo 565 !enddo 566 !enddo 567 568 call xmpi_barrier(spaceComm) 569 !if(mpi_enreg%me.eq.0) write(177,*) "end",psichi 570 !if(mpi_enreg%me.eq.1) write(178,*) "end",psichi 571 !if(mpi_enreg%me.eq.2) write(179,*) "end",psichi 572 call xmpi_barrier(spaceComm) 573 !========================================================================== 574 !********* WRITE psichi in file for reference 575 !========================================================================== 576 if(me.eq.0) then 577 call psichi_print(dtset,cryst_struc%nattyp,cryst_struc%ntypat,nkpt,my_nspinor,& 578 & nsppol,paw_dmft,pawtab,psps,t2g,x2my2d) 579 end if ! proc=me 580 581 !********************* Check normalization and occupations *************** 582 ! Only if the complete BZ is sampled (ie paw_dmft%kspectralfunc=0) 583 !========================================================================== 584 if(paw_dmft%dmft_kspectralfunc==0) then 585 ABI_MALLOC(xocc_check,(natom)) 586 ABI_MALLOC(xnorm_check,(natom)) 587 call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,xocc_check) 588 call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,xnorm_check) 589 call psichi_check(dtset,cryst_struc%nattyp,nkpt,my_nspinor,& 590 & nsppol,cryst_struc%ntypat,paw_dmft,pawtab,psps,xocc_check,xnorm_check) 591 !========================================================================== 592 !*************** write checks ******************************************* 593 !========================================================================== 594 if(abs(dtset%pawprtvol)>=3) then 595 write(message,*) "normalization computed" 596 call wrtout(std_out, message,'COLL') 597 end if 598 599 ABI_MALLOC(loc_occ_check,(natom)) 600 ABI_MALLOC(loc_norm_check,(natom)) 601 call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,loc_occ_check) 602 call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,loc_norm_check) 603 call copy_matlu(xocc_check,loc_occ_check,natom) 604 call copy_matlu(xnorm_check,loc_norm_check,natom) 605 606 write(message,'(2a,i4)') ch10," == Check: Occupations and Norm from psichi are" 607 call wrtout(std_out, message,'COLL') 608 609 if(paw_dmft%dmftcheck>=1) then 610 ! print occupations 611 write(message,'(2a,i4)') ch10,' ------ Unsymetrised Occupation' 612 call wrtout(std_out, message,'COLL') 613 614 call print_matlu(xocc_check,natom,dtset%pawprtvol) 615 616 ! print norms 617 write(message,'(2a,i4)') ch10,' ------ Unsymetrised Norm' 618 call wrtout(std_out, message,'COLL') 619 620 call print_matlu(xnorm_check,natom,dtset%pawprtvol) 621 end if 622 623 ! symetrise and print occupations 624 call sym_matlu(cryst_struc,loc_occ_check,pawang,paw_dmft) 625 626 write(message,'(2a,i4)') ch10,' ------ Symetrised Occupation' 627 call wrtout(std_out, message,'COLL') 628 629 call print_matlu(loc_occ_check,natom,dtset%pawprtvol) 630 631 ! symetrise and print norms 632 call sym_matlu(cryst_struc,loc_norm_check,pawang,paw_dmft) 633 634 write(message,'(2a,i4)') ch10,' ------ Symetrised Norm' 635 call wrtout(std_out, message,'COLL') 636 637 call print_matlu(loc_norm_check,natom,dtset%pawprtvol) 638 639 ! deallocations 640 do iatom=1,natom 641 dft_occup%matlu(iatom)%mat=loc_occ_check(iatom)%mat 642 end do 643 644 ! Tests density matrix DFT+U and density matrix computed here. 645 if(paw_dmft%dmftcheck==2.or.(paw_dmft%dmftbandi==1)) then 646 ABI_MALLOC(matlu_temp,(natom)) 647 call init_matlu(natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,matlu_temp) 648 do iatom=1,natom 649 if(paw_dmft%lpawu(iatom).ne.-1) then 650 ldim=2*paw_dmft%lpawu(iatom)+1 651 nsploop=max(paw_dmft%nsppol,paw_dmft%nspinor**2) 652 do idijeff=1,nsploop 653 ispinor=0 654 ispinor1=0 655 if(nsploop==2) then 656 isppol=spinor_idxs(1,idijeff) 657 ispinor=1 658 ispinor1=1 659 else if(nsploop==4) then 660 isppol=1 661 ispinor=spinor_idxs(1,idijeff) 662 ispinor1=spinor_idxs(2,idijeff) 663 else if(nsploop==1) then 664 isppol=1 665 ispinor=1 666 ispinor1=1 667 else 668 write(message,'(2a)') " BUG in datafordmft: nsploop should be equal to 2 or 4" 669 call wrtout(std_out,message,'COLL') 670 end if 671 do im1 = 1 , ldim 672 do im = 1 , ldim 673 if(my_nspinor==2) matlu_temp(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=& 674 & cmplx(paw_ij(iatom)%noccmmp(1,im,im1,idijeff),paw_ij(iatom)%noccmmp(2,im,im1,idijeff),kind=dp) 675 if(my_nspinor==1) matlu_temp(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=& 676 & cmplx(paw_ij(iatom)%noccmmp(1,im,im1,idijeff),zero,kind=dp) 677 end do 678 end do 679 end do 680 end if 681 end do 682 if(paw_dmft%dmftcheck==2) option=1 683 if(paw_dmft%dmftcheck<=1) option=0 684 call diff_matlu("DFT+U density matrix from INPUT wfk",& 685 & "Direct calculation of density matrix with psichi from DIAGONALIZED wfk",& 686 & matlu_temp,dft_occup%matlu,natom,option,tol3,ierrr) !tol1 tol2 tol3 687 if(ierrr==-1) then 688 write(message,'(10a)') ch10,& 689 & ' -> These two quantities should agree if three conditions are fulfilled',ch10,& 690 & ' - input wavefunctions come from the same Hamiltonien (e.g LDA/GGA)',ch10,& 691 & ' - dmatpuopt is equal to 1',ch10,& 692 & ' - all valence states are in the valence',ch10,& 693 & ' (for experts users: it is not compulsary that these conditions are fulfilled)' 694 call wrtout(std_out,message,'COLL') 695 end if 696 ! write(message,'(2a)') ch10,& 697 ! & ' ***** => Calculations of density matrices with projections and in DFT+U are coherent****' 698 ! call wrtout(std_out,message,'COLL') 699 700 call destroy_matlu(matlu_temp,natom) 701 ABI_FREE(matlu_temp) 702 else 703 write(message,'(2a)') ch10,& 704 & ' Warning: Consistency of density matrices computed from projection has not been checked: use dmftcheck>=2 ' 705 call wrtout(std_out,message,'COLL') 706 end if 707 708 call destroy_matlu(loc_norm_check,natom) 709 ABI_FREE(loc_norm_check) 710 call destroy_matlu(loc_occ_check,natom) 711 ABI_FREE(loc_occ_check) 712 713 call destroy_matlu(xocc_check,natom) 714 call destroy_matlu(xnorm_check,natom) 715 ABI_FREE(xocc_check) 716 ABI_FREE(xnorm_check) 717 endif 718 719 if(present(nbandkss)) then 720 if((me.eq.0.and.nbandkss/=0).or.(paw_dmft%dmft_kspectralfunc==1)) then 721 ! opt_renorm=1 ! if ucrpa==1, no need for individual orthonormalization 722 opt_renorm=3 723 if(dtset%ucrpa>=1.or.paw_dmft%dmft_kspectralfunc==1) opt_renorm=2 724 call psichi_renormalization(cryst_struc,paw_dmft,pawang,opt=opt_renorm) 725 call psichi_print(dtset,cryst_struc%nattyp,cryst_struc%ntypat,nkpt,my_nspinor,& 726 & nsppol,paw_dmft,pawtab,psps,t2g,x2my2d) 727 end if ! proc=me 728 end if
ABINIT/m_datafordmft [ Modules ]
NAME
m_datafordmft
FUNCTION
This module produces inputs for the DMFT calculation
COPYRIGHT
Copyright (C) 2006-2024 ABINIT group (BAmadon) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 MODULE m_datafordmft 27 28 use defs_basis 29 use defs_abitypes 30 use defs_wvltypes 31 use m_abicore 32 use m_errors 33 use m_xmpi 34 use m_dtset 35 36 use defs_datatypes, only : pseudopotential_type 37 use m_io_tools, only : open_file 38 use m_crystal, only : crystal_t 39 use m_matlu, only: matlu_type,init_matlu,sym_matlu,copy_matlu,print_matlu,diff_matlu,destroy_matlu, checkdiag_matlu, & 40 gather_matlu,add_matlu 41 use m_oper, only : init_oper,oper_type,identity_oper,loc_oper,destroy_oper,diff_oper, upfold_oper,copy_oper,prod_oper 42 use m_pawang, only : pawang_type 43 use m_pawtab, only : pawtab_type 44 use m_paw_ij, only : paw_ij_type 45 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free 46 use m_paw_dmft, only: paw_dmft_type 47 use m_mpinfo, only : proc_distrb_cycle 48 use m_matrix, only : invsqrt_matrix 49 50 implicit none 51 52 private 53 54 public :: datafordmft 55 public :: compute_levels 56 public :: psichi_renormalization 57 public :: hybridization_asymptotic_coefficient
m_datafordmft/compute_levels [ Functions ]
[ Top ] [ m_datafordmft ] [ Functions ]
NAME
compute_levels
FUNCTION
Compute levels for ctqmc
INPUTS
OUTPUT
NOTES
SOURCE
1032 subroutine compute_levels(cryst_struc,energy_level,hdc,pawang,paw_dmft,nondiag) 1033 1034 !Arguments ------------------------------------ 1035 !scalars 1036 type(crystal_t),intent(in) :: cryst_struc 1037 type(pawang_type), intent(in) :: pawang 1038 type(oper_type), intent(in) :: hdc 1039 type(paw_dmft_type), intent(in) :: paw_dmft 1040 type(oper_type),intent(inout) :: energy_level !vz_i 1041 logical, optional, intent(out) :: nondiag 1042 1043 !Local variables ------------------------------ 1044 ! scalars 1045 integer :: iatom,iband,ikpt,im1,isppol,ispinor 1046 integer :: lpawu,mbandc,natom,nspinor,nsppol,nkpt 1047 character(len=500) :: message 1048 ! arrays 1049 !************************************************************************ 1050 1051 mbandc=paw_dmft%mbandc 1052 nkpt=paw_dmft%nkpt 1053 nsppol=paw_dmft%nsppol 1054 nspinor=paw_dmft%nspinor 1055 natom=paw_dmft%natom 1056 if(present(nondiag)) nondiag=.false. 1057 1058 !======================== 1059 !Get KS eigenvalues 1060 !======================== 1061 do iband=1,mbandc 1062 do ikpt=1,nkpt 1063 do isppol=1,nsppol 1064 ! Take \epsilon_{nks} 1065 ! ======================== 1066 energy_level%ks(isppol,ikpt,iband,iband)=paw_dmft%eigen_dft(isppol,ikpt,iband) 1067 end do 1068 end do 1069 end do 1070 1071 1072 !====================================================================== 1073 !Compute atomic levels from projection of \epsilon_{nks} and symetrize 1074 !====================================================================== 1075 call loc_oper(energy_level,paw_dmft,1) 1076 ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels before sym and only DFT" 1077 ! call wrtout(std_out,message,'COLL') 1078 ! call print_matlu(energy_level%matlu,natom,1) 1079 do iatom = 1 , natom 1080 lpawu=paw_dmft%lpawu(iatom) 1081 if(lpawu/=-1) then 1082 do isppol=1,nsppol 1083 do ispinor=1,nspinor 1084 do im1=1,2*lpawu+1 1085 energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)= & 1086 & energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)& 1087 & -hdc%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)-paw_dmft%fermie 1088 end do 1089 end do 1090 end do 1091 ! write(std_out,*) "DC,fermie",hdc%matlu(iatom)%mat(1,1,1,1,1),paw_dmft%fermie 1092 end if 1093 end do ! natom 1094 call sym_matlu(cryst_struc,energy_level%matlu,pawang,paw_dmft) 1095 if(present(nondiag)) call checkdiag_matlu(energy_level%matlu,natom,tol7,nondiag) 1096 1097 write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels for Fermi Level=",paw_dmft%fermie 1098 call wrtout(std_out,message,'COLL') 1099 !call print_oper(energy_level,1,paw_dmft,1) 1100 call print_matlu(energy_level%matlu,natom,1) 1101 1102 1103 1104 end subroutine compute_levels
m_datafordmft/hybridization_asymptotic_coefficient [ Functions ]
[ Top ] [ m_datafordmft ] [ Functions ]
NAME
hybridization_asymptotic_coefficient
FUNCTION
Compute some components for the limit of hybridization
INPUTS
cryst_struc <type(crystal_t)>=crystal structure data dft_occup paw_dmft = data for self-consistent DFT+DMFT calculations. pawang <type(pawang)>=paw angular mesh and related data pawtab <type(pawtab)>
OUTPUT
paw_dmft = data for self-consistent DFT+DMFT calculations.
NOTES
SOURCE
1698 subroutine hybridization_asymptotic_coefficient(cryst_struc,paw_dmft,pawang,hybri_coeff) 1699 1700 !Arguments ------------------------------------ 1701 !scalars 1702 type(crystal_t),intent(in) :: cryst_struc 1703 type(paw_dmft_type), intent(in) :: paw_dmft 1704 type(pawang_type), intent(in) :: pawang 1705 type(matlu_type), intent(inout) :: hybri_coeff(paw_dmft%natom) 1706 !Local variables ------------------------------ 1707 type(oper_type) :: ham_a 1708 type(oper_type) :: ham_b 1709 type(oper_type) :: ham_squarelocal 1710 type(oper_type) :: ham_squareks 1711 integer :: iband1,iband2,ikpt,isppol 1712 !************************************************************************ 1713 1714 ! call init_oper(paw_dmft,self_minus_hdc) 1715 call init_oper(paw_dmft,ham_a) 1716 call init_oper(paw_dmft,ham_b) 1717 call init_oper(paw_dmft,ham_squareks) 1718 call init_oper(paw_dmft,ham_squarelocal) 1719 1720 ! Create self_minus_hdc%matlu = Sigma-Hdc in local orbitals 1721 ! call add_matlu(self%oper(paw_dmft%dmft_nwlo)%matlu,self%hdc%matlu,& 1722 !& self_minus_hdc%matlu,cryst_struc%natom,-1) 1723 1724 ! write(message,'(a,2x,a)') ch10, " == self_minus_hdc (1)" 1725 ! call wrtout(std_out,message,'COLL') 1726 ! call print_matlu(self_minus_hdc%matlu,paw_dmft%natom,1,opt_exp=1) 1727 1728 !! Create self_minus_hdc%ks = Upfold Sigma-Hdc 1729 ! call upfold_oper(self_minus_hdc,paw_dmft,1) 1730 ! call loc_oper(self_minus_hdc,paw_dmft,1) 1731 1732 ! write(message,'(a,2x,a)') ch10, " == self_minus_hdc (2)" 1733 ! call wrtout(std_out,message,'COLL') 1734 ! call print_matlu(self_minus_hdc%matlu,paw_dmft%natom,1,opt_exp=1) 1735 1736 ! Create ham_a%ks = H_ks in KS basis 1737 !---------------------------------------------------- 1738 do iband1=1,paw_dmft%mbandc 1739 do iband2=1,paw_dmft%mbandc 1740 do ikpt=1,paw_dmft%nkpt 1741 do isppol=1,paw_dmft%nsppol 1742 if(iband1==iband2) then 1743 ham_a%ks(isppol,ikpt,iband1,iband2) = cmplx(paw_dmft%eigen_dft(isppol,ikpt,iband1),0.d0,kind=dp) 1744 else 1745 ham_a%ks(isppol,ikpt,iband1,iband2) = czero 1746 end if 1747 end do 1748 end do 1749 end do 1750 end do 1751 1752 ! Create ham_a%matlu = H_ks in local orbitals 1753 !--------------------------------------------- 1754 call loc_oper(ham_a,paw_dmft,1) 1755 1756 ! Symetrise the local quantity (energy levels) 1757 !--------------------------------------------- 1758 call sym_matlu(cryst_struc,ham_a%matlu,pawang,paw_dmft) 1759 1760 ! Create ham_b%ks : Duplicate both ks and local part of ham_a into ham_b 1761 !----------------------------------------------------------------------- 1762 call copy_oper(ham_a,ham_b) 1763 1764 ! Compute ham_squareks%ks : Make a product of the two KS version 1765 !------------------------------------------------------------------ 1766 call prod_oper(ham_a,ham_b,ham_squareks,1) 1767 1768 ! Compute ham_squareks%matlu 1769 !--------------------------- 1770 call loc_oper(ham_squareks,paw_dmft,1) 1771 1772 ! Symetrise ham_squareks%matlu 1773 !------------------------------ 1774 call sym_matlu(cryst_struc,ham_squareks%matlu,pawang,paw_dmft) 1775 1776 ! write(message,'(a,2x,a)') ch10, " == squareks" 1777 ! call wrtout(std_out,message,'COLL') 1778 ! call print_matlu(ham_squareks%matlu,paw_dmft%natom,1,opt_exp=1) 1779 1780 1781 ! Compute ham_squarelocal%matlu 1782 !------------------------------- 1783 call prod_oper(ham_a,ham_b,ham_squarelocal,2) 1784 1785 ! Compute the product in local orbitals 1786 !-------------------------------------- 1787 call sym_matlu(cryst_struc,ham_squarelocal%matlu,pawang,paw_dmft) 1788 1789 ! write(message,'(a,2x,a)') ch10, " == squarelocal" 1790 ! call wrtout(std_out,message,'COLL') 1791 ! call print_matlu(ham_squarelocal%matlu,paw_dmft%natom,1,opt_exp=1) 1792 1793 ! The difference of ham_squareks and ham_squarelocal 1794 ! gives the coefficient that we are looking for ( such that F_ij(iw_n) = -C_ij/(iw_n) ). 1795 !---------------------------------------------------------------------------------------- 1796 call add_matlu(ham_squareks%matlu,ham_squarelocal%matlu,hybri_coeff,cryst_struc%natom,-1) 1797 1798 ! write(message,'(a,2x,a)') ch10, " == Coeff C_ij before sym" 1799 ! call wrtout(std_out,message,'COLL') 1800 ! call print_matlu(hybri_coeff,paw_dmft%natom,1,opt_exp=1) 1801 1802 ! Symetrise the local quantity 1803 !------------------------------ 1804 call sym_matlu(cryst_struc,hybri_coeff,pawang,paw_dmft) 1805 1806 ! write(message,'(a,2x,a)') ch10, " == Coeff C_ij after sym" 1807 ! call wrtout(std_out,message,'COLL') 1808 ! call print_matlu(hybri_coeff,paw_dmft%natom,1,opt_exp=1) 1809 1810 call destroy_oper(ham_squarelocal) 1811 ! call destroy_oper(self_minus_hdc) 1812 call destroy_oper(ham_a) 1813 call destroy_oper(ham_b) 1814 call destroy_oper(ham_squareks) 1815 1816 1817 end subroutine hybridization_asymptotic_coefficient
m_datafordmft/psichi_check [ Functions ]
[ Top ] [ m_datafordmft ] [ Functions ]
NAME
psichi_check
FUNCTION
Check psichi: compute occupations
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset maxnproju = maximum number of projector for DFT+U species nattyp(ntypat)= # atoms of each type mband= number of bands nkpt=number of k points my_nspinor=number of spinorial components of the wavefunctions (on current proc) nsppol=1 for unpolarized, 2 for spin-polarized ntypat= number of species paw_dmft <type(paw_dmft)>=paw data for the self-consistency pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials OUTPUTS: xocc_check(nsppol,my_nspinor,my_nspinor,natom,2*maxlpawu+1,2*maxlpawu+1): density matrix xnorm_check(nsppol,my_nspinor,my_nspinor,natom,2*maxlpawu+1,2*maxlpawu+1): matrix of norms
SIDE EFFECTS
check psichi: compute norm and occupations
SOURCE
922 subroutine psichi_check(dtset,nattyp,nkpt,my_nspinor,& 923 & nsppol,ntypat,paw_dmft,pawtab,psps,xocc_check,xnorm_check) 924 925 !Arguments ------------------------------------ 926 !scalars 927 integer,intent(in) :: nkpt,my_nspinor,nsppol,ntypat 928 !arrays 929 integer, intent(in) :: nattyp(ntypat) 930 type(dataset_type),intent(in) :: dtset 931 type(pseudopotential_type),intent(in) :: psps 932 type(paw_dmft_type), intent(in) :: paw_dmft 933 type(matlu_type), intent(inout) :: xocc_check(dtset%natom) !vz_i 934 type(matlu_type), intent(inout) :: xnorm_check(dtset%natom) !vz_i 935 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 936 !Local variables ------------------------------------ 937 integer :: band_index,facnsppol,ibg,isppol,ikpt,iband,ibandc,ispinor,icat,itypat 938 integer :: iat,iatom,ilmn,lmn_size,m,m1,nband_k 939 complex(dpc) :: psichic,psichic1 940 real(dp) :: chinorm 941 942 ! ********************************************************************* 943 facnsppol=1 944 if(my_nspinor==1.and.nsppol==1) then 945 facnsppol=2 946 end if 947 948 ibg=0 949 band_index=0 950 do isppol=1,nsppol 951 do ikpt=1,nkpt 952 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 953 ibandc=0 954 do iband=1,nband_k 955 if(paw_dmft%band_in(iband)) ibandc=ibandc+1 956 do ispinor=1,my_nspinor 957 icat=1 958 iat=0 959 do itypat=1,dtset%ntypat 960 lmn_size=pawtab(itypat)%lmn_size 961 do iatom=icat,icat+nattyp(itypat)-1 962 iat=iat+1 963 ! ------------ Select correlated atoms 964 if(paw_dmft%lpawu(iatom).ne.-1) then 965 chinorm=1.d0 966 do ilmn=1,lmn_size 967 ! ------------ Select l=lpawu. 968 if (psps%indlmn(1,ilmn,itypat)==paw_dmft%lpawu(iatom).and.& 969 & psps%indlmn(3,ilmn,itypat)==1) then 970 do ilmn1=1,lmn_size 971 ! ------------ Select l=lpawu and do not sum over projectors 972 ! (this is already done in paw_dmft%psichi) 973 if (psps%indlmn(1,ilmn1,itypat)==paw_dmft%lpawu(iatom).and.& 974 & psps%indlmn(3,ilmn1,itypat)==1) then 975 ! ------------ Check that the band is choosen (within nbandi and nbandf) 976 if(paw_dmft%band_in(iband)) then 977 m=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1 978 m1=psps%indlmn(2,ilmn1,itypat)+psps%indlmn(1,ilmn,itypat)+1 979 if(psps%indlmn(3,ilmn,itypat)==1) then 980 do ispinor1=1,my_nspinor 981 psichic=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m) 982 psichic1=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor1,iat,m1) 983 ! ------------ Compute occupation matrix 984 xocc_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)=& 985 & xocc_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)& 986 ! & +conjg(psichic)*psichic1*dtset%wtk(ikpt)*facpara*occ(iband+band_index) 987 & +conjg(psichic1)*psichic*dtset%wtk(ikpt)*facpara*occ(iband+band_index)/facnsppol 988 ! ------------ Compute norm (should be equal to noccmmp 989 ! (dmatpuopt=1) if all bands are taken into account) 990 xnorm_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)=& 991 & xnorm_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)& 992 ! & +conjg(psichic)*psichic1*dtset%wtk(ikpt)*facpara 993 & +conjg(psichic1)*psichic*dtset%wtk(ikpt)*facpara 994 end do ! ispinor1 995 end if 996 end if ! paw_dmft%band_in 997 end if 998 end do !ilmn1 999 end if 1000 end do !ilmn 1001 end if ! lpawu.ne.-1 1002 end do ! iatom 1003 icat=icat+nattyp(itypat) 1004 end do ! itypat 1005 end do ! ispinor 1006 end do !iband 1007 band_index=band_index+nband_k 1008 ibg=ibg+nband_k*my_nspinor 1009 end do !ikpt 1010 end do ! isppol 1011 1012 end subroutine psichi_check 1013 !DBG_EXIT("COLL") 1014 end subroutine datafordmft
m_datafordmft/psichi_print [ Functions ]
[ Top ] [ m_datafordmft ] [ Functions ]
NAME
psichi_print
FUNCTION
Print psichi for reference
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset maxnproju = maximum number of projector for DFT+U species nattyp(ntypat)= # atoms of each type mband= number of bands nkpt=number of k points my_nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized paw_dmft <type(paw_dmft)>=paw data for the self-consistency pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psichi(2,nsppol,nkpt,mband,my_nspinor,dtset%natom,(2*paw_dmft%maxlpawu+1))) projections for DMFT psps <type(pseudopotential_type)>=variables related to pseudopotentials
SIDE EFFECTS
print psichi in forlb.ovlp
SOURCE
759 subroutine psichi_print(dtset,nattyp,ntypat,nkpt,my_nspinor,& 760 &nsppol,paw_dmft,pawtab,psps,t2g,x2my2d) 761 762 !Arguments ------------------------------------ 763 !scalars 764 integer,intent(in) :: nkpt,my_nspinor,nsppol,ntypat 765 !arrays 766 integer, intent(in) :: nattyp(ntypat) 767 type(dataset_type),intent(in) :: dtset 768 type(pseudopotential_type),intent(in) :: psps 769 logical t2g,x2my2d 770 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 771 type(paw_dmft_type), intent(in) :: paw_dmft 772 !Local variables ------------------------------------ 773 integer :: ibg,isppol,ikpt,iband,ibandc,ispinor,icat,itypat,lmn_size 774 integer :: iat,iatom,jj1,ilmn,m1,nband_k,unt 775 integer :: m1_t2g,ll,m1_x2my2d 776 real(dp) :: chinorm 777 character(len=500) :: msg 778 779 ! ********************************************************************* 780 ll=1 781 782 if (open_file('forlb.ovlp',msg,newunit=unt,form='formatted',status='unknown') /= 0) then 783 ABI_ERROR(msg) 784 end if 785 rewind(unt) 786 787 ! Header for calc_uCRPA.F90 788 write(unt,*) "# isppol nspinor natom m Re(<psi|chi>) Im(<psi|chi>)" 789 if (COUNT(pawtab(:)%lpawu.NE.-1).EQ.1) then 790 do itypat=1,ntypat 791 if(t2g) then 792 if(pawtab(itypat)%lpawu.ne.-1) write(unt,*) "l= ",ll,itypat 793 else if(x2my2d) then 794 if(pawtab(itypat)%lpawu.ne.-1) write(unt,*) "l= ",ll-1,itypat 795 else 796 if(pawtab(itypat)%lpawu.ne.-1) write(unt,*) "l= ",pawtab(itypat)%lpawu,itypat 797 end if 798 end do 799 else 800 write(unt,*) "More than one correlated species" 801 end if 802 803 write(unt,*) "Bands ",paw_dmft%dmftbandi,paw_dmft%dmftbandf 804 805 ibg=0 806 do isppol=1,nsppol 807 do ikpt=1,nkpt 808 ! rewind(1023) 809 write(unt,'(a6,2x,i6)') "ikpt =",ikpt 810 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 811 ibandc=0 812 do iband=1,nband_k 813 if(paw_dmft%band_in(iband)) then 814 ibandc=ibandc+1 815 write(unt,'(a8,2x,i6)') " iband =",iband 816 end if 817 do ispinor=1,my_nspinor 818 icat=1 819 ! write(std_out,*) isppol,ikpt,iband,ispinor 820 iat=0 ! to declare 821 do itypat=1,dtset%ntypat 822 lmn_size=pawtab(itypat)%lmn_size 823 ! write(std_out,*) isppol,ikpt,iband,ispinor 824 do iatom=icat,icat+nattyp(itypat)-1 825 iat=iat+1 826 jj1=0 827 if(pawtab(itypat)%lpawu.ne.-1) then 828 ! chinorm=(pawtab(itypat)%phiphjint(1)) 829 chinorm=1.d0 830 ! write(std_out,*) isppol,ikpt,iband,ispinor,iat 831 m1_t2g=0 832 m1_x2my2d=0 833 do ilmn=1,lmn_size 834 ! write(std_out,*) ilmn 835 ! ------------ Select l=lpawu. --------------------------------------- 836 if (psps%indlmn(1,ilmn,itypat)==pawtab(itypat)%lpawu.and.psps%indlmn(3,ilmn,itypat)==1) then 837 ! ------------ Check that the band is choosen (within nbandi and nbandf) 838 if(paw_dmft%band_in(iband)) then 839 jj1=jj1+1 840 if(jj1>(2*pawtab(itypat)%lpawu+1)) then 841 write(message,'(a,a,i4,i5,i4)') ch10," Error 2 in datafordmft",jj1,pawtab(itypat)%lpawu 842 call wrtout(std_out, message,'COLL') 843 ABI_ERROR("Aborting now") 844 end if ! jj1 845 ! if(jj1>pawtab(dtset%typat(iatom))%nproju*(2*lpawu+1)) then 846 ! write(message,'(a,a,a,a)') ch10,& 847 ! & 'BUG: jj1 is not correct in datafordmft psichi_print',ch10,& 848 ! & 'Action: CONTACT Abinit group' 849 ! call wrtout(std_out, message,'COLL') 850 ! call abi_abort('COLL') 851 ! end if ! jj1 852 m1=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1 853 ! ----- Print only when the sum over projectors is done 854 ! write(std_out,*) ilmn,m1 855 if(t2g) then 856 if(m1==1.or.m1==2.or.m1==4) then 857 m1_t2g=m1_t2g+1 858 write(unt,'(4i6,3x,2f23.15)') isppol, ispinor, iat, m1,& 859 & real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g))/chinorm,& 860 & aimag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g))/chinorm 861 end if 862 else if(x2my2d) then 863 if(m1==5) then 864 m1_x2my2d=1 865 write(unt,'(4i6,3x,2f23.15)') isppol, ispinor, iat, m1,& 866 & real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_x2my2d))/chinorm,& 867 & aimag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_x2my2d))/chinorm 868 end if 869 else 870 write(unt,'(4i6,3x,2f23.15)') isppol, ispinor, iat, m1,& 871 & real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1))/chinorm,& 872 & aimag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1))/chinorm 873 end if 874 ! if natom=1 then jj1 maximal value should be 2*lpawu+1 875 end if ! paw_dmft%band_in 876 end if 877 end do !ilmn 878 end if ! lpawu.ne.-1 879 end do ! iatom 880 icat=icat+nattyp(itypat) 881 end do ! itypat 882 end do ! ispinor 883 end do !iband 884 ibg=ibg+nband_k*my_nspinor 885 end do !ikpt 886 end do ! isppol 887 ! write(unt,*) "Fermi level (in Ryd)=" 888 ! write(unt,*) fermie*two 889 close(unt) 890 end subroutine psichi_print
m_datafordmft/psichi_renormalization [ Functions ]
[ Top ] [ m_datafordmft ] [ Functions ]
NAME
psichi_renormalization
FUNCTION
Renormalize psichi.
INPUTS
cryst_struc <type(crystal_t)>= crystal structure data. paw_dmft = data for DFT+DMFT calculations. pawang <type(pawang)>=paw angular mesh and related data
OUTPUT
paw_dmft%psichi(nsppol,nkpt,mband,nspinor,dtset%natom,(2*maxlpawu+1))): projections <Psi|chi> are orthonormalized.
NOTES
SOURCE
1127 subroutine psichi_renormalization(cryst_struc,paw_dmft,pawang,opt) 1128 1129 !Arguments ------------------------------------ 1130 !scalars 1131 type(crystal_t),intent(in) :: cryst_struc 1132 type(paw_dmft_type), intent(inout) :: paw_dmft 1133 type(pawang_type), intent(in) :: pawang 1134 integer, optional, intent(in) :: opt 1135 1136 !Local variables ------------------------------ 1137 !scalars 1138 integer :: iatom,ib,ikpt,im,ispinor,isppol,jkpt 1139 integer :: natom,mbandc,ndim,nkpt,nspinor,nsppol,option 1140 real(dp), pointer :: temp_wtk(:) => null() 1141 real(dp) :: pawprtvol 1142 type(oper_type) :: norm 1143 type(oper_type) :: oper_temp 1144 character(len=500) :: message 1145 !arrays 1146 ! real(dp),allocatable :: e0pde(:,:,:),omegame0i(:) 1147 !************************************************************************ 1148 1149 DBG_ENTER("COLL") 1150 1151 option=3 1152 if(present(opt)) then 1153 if(opt==2) option=2 1154 if(opt==3) option=3 1155 end if 1156 pawprtvol=2 1157 1158 nsppol = paw_dmft%nsppol 1159 nkpt = paw_dmft%nkpt 1160 mbandc = paw_dmft%mbandc 1161 natom = cryst_struc%natom 1162 nspinor = paw_dmft%nspinor 1163 1164 1165 !== Normalize psichi 1166 if(option==1) then 1167 ! ==================================== 1168 ! == simply renormalize psichi ======= 1169 ! ==================================== 1170 write(message,'(2a)') ch10," Psichi are renormalized " 1171 call wrtout(std_out, message,'COLL') 1172 do isppol=1,nsppol 1173 do ikpt=1,nkpt 1174 do ib=1,mbandc 1175 do iatom=1,natom 1176 if(paw_dmft%lpawu(iatom).ne.-1) then 1177 ndim=2*paw_dmft%lpawu(iatom)+1 1178 do im=1,ndim 1179 do ispinor=1,nspinor 1180 ! write(std_out,*) "psichi1",paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im) 1181 paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im)= & 1182 & paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im)/ & 1183 & sqrt(real(norm%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor))) 1184 end do ! ispinor 1185 end do ! im 1186 end if 1187 end do ! iatom 1188 end do ! ib 1189 end do ! ikpt 1190 end do ! isppol 1191 ! todo_ab introduce correct orthonormalization in the general case. 1192 1193 else if(option==2) then ! option==2 1194 ! ==================================== 1195 ! == renormalize k-point after k-point 1196 ! ==================================== 1197 ABI_MALLOC(temp_wtk,(1)) 1198 1199 write(message,'(2a,i5)') ch10,' Nb of K-point',nkpt 1200 call wrtout(std_out,message,'COLL') 1201 do jkpt=1,nkpt ! jkpt 1202 write(message,'(2a,i5)') ch10,' == Renormalization for K-point',jkpt 1203 call wrtout(std_out,message,'COLL') 1204 temp_wtk(1)=one 1205 1206 call normalizepsichi(cryst_struc,1,paw_dmft,pawang,temp_wtk,jkpt) 1207 end do ! jkpt 1208 write(message,'(2a)') ch10,' ===== K-points all renormalized' 1209 call wrtout(std_out,message,'COLL') 1210 ABI_FREE(temp_wtk) 1211 1212 else if(option==3) then ! option==3 1213 ! ==================================== 1214 ! == renormalize the sum over k-points 1215 ! ==================================== 1216 ! allocate(temp_wtk(nkpt)) 1217 temp_wtk=>paw_dmft%wtk 1218 write(message,'(6a)') ch10,' ====================================== '& 1219 & ,ch10,' == Renormalization for all K-points == '& 1220 & ,ch10,' =======================================' 1221 call wrtout(std_out,message,'COLL') 1222 call normalizepsichi(cryst_struc,nkpt,paw_dmft,pawang,temp_wtk) 1223 ! deallocate(temp_wtk) 1224 1225 end if 1226 1227 !== Change back repr for norm 1228 1229 1230 !=============================================== 1231 !== Compute norm with new psichi 1232 !=============================================== 1233 write(message,'(2a)') ch10,' ===== Compute norm with new psichi' 1234 call wrtout(std_out,message,'COLL') 1235 call init_oper(paw_dmft,norm,nkpt=paw_dmft%nkpt,wtk=paw_dmft%wtk) 1236 !== Build identity for norm%ks (option=1) 1237 call identity_oper(norm,1) 1238 1239 !== Deduce norm%matlu from norm%ks with new psichi 1240 if (paw_dmft%dmft_kspectralfunc==1) then 1241 call init_oper(paw_dmft,oper_temp) 1242 do jkpt=1,nkpt ! jkpt 1243 call loc_oper(norm,paw_dmft,1,jkpt=jkpt) 1244 write(message,'(2a,i5)') & 1245 & ch10," == Check: Overlap with renormalized psichi for k-point",jkpt 1246 call wrtout(std_out,message,'COLL') 1247 call print_matlu(norm%matlu,natom,prtopt=1) 1248 !== Check that norm is now the identity 1249 call identity_oper(oper_temp,2) 1250 call diff_matlu('Overlap after renormalization','Identity',& 1251 & norm%matlu,oper_temp%matlu,norm%natom,1,tol6,zero_or_one=1) 1252 enddo 1253 call destroy_oper(oper_temp) 1254 else !dmft_kspectralfunc 1255 write(message,'(2a)') ch10,' ===== Calling loc_oper' 1256 call wrtout(std_out,message,'COLL') 1257 call loc_oper(norm,paw_dmft,option) 1258 write(message,'(2a)') ch10,' ===== Finished loc_oper' 1259 call wrtout(std_out,message,'COLL') 1260 1261 1262 !== Print norm%matlu unsymetrized with new psichi 1263 if(pawprtvol>2) then 1264 write(message,'(4a,2a)') & 1265 & ch10," == Check: Overlap with renormalized psichi without symetrization is == " 1266 call wrtout(std_out,message,'COLL') 1267 call print_matlu(norm%matlu,natom,prtopt=1) 1268 end if 1269 1270 1271 !== Symetrise norm%matlu with new psichi 1272 call sym_matlu(cryst_struc,norm%matlu,pawang,paw_dmft) 1273 1274 !== Print norm%matlu symetrized with new psichi 1275 if(pawprtvol>2) then 1276 write(message,'(4a,2a)') & 1277 & ch10," == Check: Overlap with renormalized psichi and symetrization is ==" 1278 call wrtout(std_out,message,'COLL') 1279 call print_matlu(norm%matlu,natom,prtopt=1,opt_diag=-1) 1280 end if 1281 1282 !== Check that norm is now the identity 1283 call init_oper(paw_dmft,oper_temp) 1284 call identity_oper(oper_temp,2) 1285 call diff_oper('Overlap after renormalization','Identity',& 1286 & norm,oper_temp,option,tol6) 1287 call destroy_oper(oper_temp) 1288 1289 call destroy_oper(norm) 1290 endif ! dmft_kspectralfunc 1291 1292 paw_dmft%lpsichiortho=1 1293 1294 DBG_EXIT("COLL") 1295 1296 CONTAINS 1297 !===========================================================
psichi_renormalization/normalizepsichi [ Functions ]
[ Top ] [ psichi_renormalization ] [ Functions ]
NAME
normalizepsichi
FUNCTION
normalizepsichi psichi from norm
INPUTS
SIDE EFFECTS
change psichi: normalizepsichi it
SOURCE
1314 subroutine normalizepsichi(cryst_struc,nkpt,paw_dmft,pawang,temp_wtk,jkpt) 1315 1316 !Arguments ------------------------------------ 1317 integer,intent(in) :: nkpt 1318 integer,optional,intent(in) :: jkpt 1319 real(dp),pointer :: temp_wtk(:) 1320 !scalars 1321 type(crystal_t),intent(in) :: cryst_struc 1322 type(paw_dmft_type), intent(inout) :: paw_dmft 1323 type(pawang_type), intent(in) :: pawang 1324 1325 !Local variables ------------------------------ 1326 !scalars 1327 integer :: diag,iatom,ib,ikpt1,im,im1,ispinor,ispinor1,isppol,isppol1,jc,jc1 1328 integer :: tndim,dum 1329 integer :: natom,mbandc,ndim,nspinor,nsppol,iortho,natomcor 1330 integer :: itot,itot1,dimoverlap,iatomcor 1331 real(dp) :: pawprtvol 1332 type(oper_type) :: norm1,norm2,norm3 1333 character(len=500) :: message 1334 complex(dpc),allocatable :: wan(:,:,:),sqrtmatinv(:,:),wanall(:) 1335 type(coeff2c_type), allocatable :: overlap(:) 1336 complex(dpc), allocatable :: largeoverlap(:,:) 1337 complex(dpc), allocatable :: psichivect(:,:) 1338 !arrays 1339 ! real(dp),allocatable :: e0pde(:,:,:),omegame0i(:) 1340 !************************************************************************ 1341 nsppol = paw_dmft%nsppol 1342 mbandc = paw_dmft%mbandc 1343 natom = cryst_struc%natom 1344 nspinor = paw_dmft%nspinor 1345 pawprtvol=3 1346 diag=0 1347 natomcor=0 1348 dimoverlap=0 1349 do iatom=1,natom 1350 if(paw_dmft%lpawu(iatom).ne.-1) then 1351 natomcor=natomcor+1 1352 ndim=2*paw_dmft%lpawu(iatom)+1 1353 tndim=nspinor*ndim 1354 dimoverlap=dimoverlap+tndim 1355 ! write(6,*) "atom, dimoverlap",iatom,dimoverlap,natomcor 1356 end if 1357 end do 1358 1359 if(nkpt/=1.and.present(jkpt)) then 1360 message = 'BUG in psichi_normalization' 1361 ABI_ERROR(message) 1362 end if 1363 1364 iortho=1 1365 ! write(6,*) "nkpt, iortho",nkpt,iortho 1366 !if (natomcor>1) iortho=2 1367 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1368 ! First case: usual case (should in fact be used only for one atom and nkpt=1) 1369 if((.not.present(jkpt)).and.iortho==1) then 1370 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1371 ! ********************************************************************* 1372 call init_oper(paw_dmft,norm1,nkpt=nkpt,wtk=temp_wtk) 1373 1374 ! == Build identity for norm1%ks (option=1) 1375 call identity_oper(norm1,1) 1376 1377 if(nkpt==1.and.present(jkpt)) then 1378 call loc_oper(norm1,paw_dmft,1,jkpt=jkpt) 1379 end if 1380 if(.not.present(jkpt)) then 1381 call loc_oper(norm1,paw_dmft,1) 1382 end if 1383 if(nkpt>1) then 1384 call sym_matlu(cryst_struc,norm1%matlu,pawang,paw_dmft) 1385 end if 1386 1387 if(pawprtvol>2) then 1388 write(message,'(2a)') ch10,' - Print norm with current psichi ' 1389 call wrtout(std_out,message,'COLL') 1390 call print_matlu(norm1%matlu,natom,prtopt=1,opt_exp=1) 1391 end if 1392 ! ==------------------------------------- 1393 ! == Start loop on atoms 1394 ABI_MALLOC(overlap,(natom)) 1395 do iatom=1,natom 1396 if(paw_dmft%lpawu(iatom).ne.-1) then 1397 ndim=2*paw_dmft%lpawu(iatom)+1 1398 tndim=nsppol*nspinor*ndim 1399 ABI_MALLOC(overlap(iatom)%value,(tndim,tndim)) 1400 overlap(iatom)%value=czero 1401 end if 1402 end do 1403 ! ==------------------------------------- 1404 1405 ! built large overlap matrix 1406 write(message,'(2a)') ch10,' - Overlap (before orthonormalization) -' 1407 call wrtout(std_out,message,'COLL') 1408 call gather_matlu(norm1%matlu,overlap,cryst_struc%natom,option=1,prtopt=1) 1409 call destroy_oper(norm1) 1410 1411 1412 1413 do iatom=1,natom 1414 if(paw_dmft%lpawu(iatom).ne.-1) then 1415 ndim=2*paw_dmft%lpawu(iatom)+1 1416 tndim=nsppol*nspinor*ndim 1417 ABI_MALLOC(sqrtmatinv,(tndim,tndim)) 1418 1419 ! == Compute Inverse Square root of overlap : O^{-0.5} 1420 !do im=1,tndim 1421 ! do im1=1,tndim 1422 ! !write(message,'(a,1x,a,e21.14,a,e21.14,a)') "overlap", & 1423 ! !"(",real(overlap(1)%value(im,im1)),",",aimag(overlap(1)%value(im,im1)),")" 1424 ! write(6,*) "overlap",overlap(iatom)%value(im,im1) 1425 ! enddo 1426 !enddo 1427 !stop 1428 !call wrtout(std_out,message,'COLL') 1429 if(diag==0) then 1430 call invsqrt_matrix(overlap(iatom)%value,tndim,dum) 1431 sqrtmatinv=overlap(iatom)%value 1432 else 1433 sqrtmatinv(:,:)=czero 1434 do ib=1,tndim 1435 sqrtmatinv(ib,ib)=cone/(sqrt(overlap(iatom)%value(ib,ib))) 1436 end do 1437 end if 1438 1439 ! == Apply O^{-0.5} on psichi 1440 ABI_MALLOC(wan,(nsppol,nspinor,ndim)) 1441 ! write(std_out,*) mbandc,nsppol,nspinor,ndim 1442 ! write(std_out,*) paw_dmft%psichi(1,1,1,1,1,1) 1443 do ikpt=1,nkpt 1444 do ib=1,mbandc 1445 if(present(jkpt)) then 1446 ikpt1=jkpt 1447 else 1448 ikpt1=ikpt 1449 end if 1450 jc=0 1451 wan=czero 1452 do isppol=1,nsppol 1453 do ispinor=1,nspinor 1454 do im=1,ndim 1455 ! write(std_out,*) "psichi", paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im) 1456 jc=jc+1 1457 jc1=0 1458 do isppol1=1,nsppol 1459 do ispinor1=1,nspinor 1460 do im1=1,ndim 1461 jc1=jc1+1 1462 wan(isppol,ispinor,im)= wan(isppol,ispinor,im) & 1463 & + paw_dmft%psichi(isppol1,ikpt1,ib,ispinor1,iatom,im1)*sqrtmatinv(jc,jc1) 1464 end do ! ispinor1 1465 end do ! isppol1 1466 end do ! im1 1467 end do ! im 1468 end do ! ispinor 1469 end do ! isppol 1470 do isppol=1,nsppol 1471 do ispinor=1,nspinor 1472 do im=1,ndim 1473 paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im)=wan(isppol,ispinor,im) 1474 ! write(std_out,*) "psichi2", paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im) 1475 end do ! ispinor 1476 end do ! isppol 1477 end do ! im 1478 end do ! ib 1479 end do ! ikpt 1480 ABI_FREE(wan) 1481 ABI_FREE(sqrtmatinv) 1482 ! write(std_out,*) paw_dmft%psichi(1,1,1,1,1,1) 1483 1484 ! ==------------------------------------- 1485 end if ! lpawu.ne.-1 1486 end do ! iatom 1487 ! == End loop on atoms 1488 ! ==------------------------------------- 1489 do iatom=1,natom 1490 if(paw_dmft%lpawu(iatom).ne.-1) then 1491 ABI_FREE(overlap(iatom)%value) 1492 end if 1493 end do 1494 ABI_FREE(overlap) 1495 1496 ! ====================================================================== 1497 ! == Check norm with new psichi. 1498 ! ====================================================================== 1499 1500 call init_oper(paw_dmft,norm1,nkpt=nkpt,wtk=temp_wtk) 1501 1502 call identity_oper(norm1,1) 1503 1504 if(nkpt==1.and.present(jkpt)) then 1505 call loc_oper(norm1,paw_dmft,1,jkpt=jkpt) 1506 end if 1507 if(.not.present(jkpt)) then 1508 call loc_oper(norm1,paw_dmft,1) 1509 end if 1510 1511 if (nkpt>1) then 1512 call sym_matlu(cryst_struc,norm1%matlu,pawang,paw_dmft) 1513 end if 1514 1515 if(pawprtvol>2) then 1516 write(message,'(2a)') ch10,' - Print norm with new psichi ' 1517 call wrtout(std_out,message,'COLL') 1518 call print_matlu(norm1%matlu,natom,prtopt=1) 1519 end if 1520 1521 ! ====================================================================== 1522 ! == Check that norm-identity is zero 1523 ! ====================================================================== 1524 call init_oper(paw_dmft,norm2,nkpt=nkpt,wtk=temp_wtk) 1525 call init_oper(paw_dmft,norm3,nkpt=nkpt,wtk=temp_wtk) 1526 call identity_oper(norm2,2) 1527 call add_matlu(norm1%matlu,norm2%matlu,norm3%matlu,natom,-1) 1528 call destroy_oper(norm2) 1529 if(pawprtvol>2) then 1530 write(message,'(2a)') ch10,' - Print norm with new psichi minus Identity ' 1531 call wrtout(std_out,message,'COLL') 1532 call print_matlu(norm3%matlu,natom,prtopt=1,opt_exp=1) 1533 end if 1534 call destroy_oper(norm3) 1535 1536 call destroy_oper(norm1) 1537 ! call flush(std_out) ! debug debug debug debug 1538 ! ABI_ERROR("Stop for debugging") 1539 1540 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1541 ! New implementation, several atoms, general case. 1542 else if(present(jkpt).or.iortho==2) then 1543 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1544 1545 ABI_MALLOC(largeoverlap,(dimoverlap,dimoverlap)) 1546 ABI_MALLOC(psichivect,(mbandc,dimoverlap)) 1547 ABI_MALLOC(sqrtmatinv,(dimoverlap,dimoverlap)) 1548 ABI_MALLOC(wanall,(dimoverlap)) 1549 1550 1551 ! Big loop over isppol 1552 do isppol=1,nsppol 1553 1554 do ib=1,mbandc 1555 itot=0 1556 do iatom=1,natom 1557 if(paw_dmft%lpawu(iatom).ne.-1) then 1558 ndim=2*paw_dmft%lpawu(iatom)+1 1559 do im=1,ndim 1560 do ispinor=1,nspinor 1561 itot=itot+1 1562 if(itot>dimoverlap) write(std_out,*) "itot>ndim",itot,ndim 1563 ! write(6,*) "ib,iatom,im,ispinor",ib,iatom,im,ispinor,jkpt 1564 psichivect(ib,itot)= paw_dmft%psichi(isppol,jkpt,ib,ispinor,iatom,im) 1565 enddo ! ispinor 1566 enddo ! im 1567 endif 1568 enddo ! iatom 1569 enddo ! ib 1570 1571 1572 ! calculation of overlap 1573 largeoverlap=czero 1574 do ib=1,mbandc 1575 do itot=1,dimoverlap 1576 do itot1=1,dimoverlap 1577 largeoverlap(itot,itot1)=largeoverlap(itot,itot1)+ & 1578 & psichivect(ib,itot)*conjg(psichivect(ib,itot1)) 1579 enddo ! itot1 1580 enddo ! itot 1581 enddo ! ib 1582 1583 ! Math: orthogonalisation of overlap 1584 write(std_out,*)"jkpt=",jkpt 1585 do itot=1,dimoverlap 1586 write(std_out,'(100f7.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap) 1587 enddo 1588 call invsqrt_matrix(largeoverlap,dimoverlap,dum) 1589 sqrtmatinv=largeoverlap 1590 write(std_out,*)"jkpt=",jkpt 1591 do itot=1,dimoverlap 1592 write(std_out,'(100f7.3)') (sqrtmatinv(itot,itot1),itot1=1,dimoverlap) 1593 enddo 1594 write(std_out,*)"jkpt=",jkpt 1595 do itot=1,dimoverlap 1596 write(std_out,'(100e9.3)') (sqrtmatinv(itot,itot1),itot1=1,dimoverlap) 1597 enddo 1598 1599 1600 do ib=1,mbandc 1601 wanall=czero 1602 do itot=1,dimoverlap 1603 do itot1=1,dimoverlap 1604 wanall(itot)= wanall(itot)+psichivect(ib,itot1)*sqrtmatinv(itot,itot1) 1605 enddo ! itot1 1606 ! write(std_out,'(3i3,2x,i3,2x,2e15.5,2x,2e15.5)') jkpt,isppol,ib,itot,psichivect(ib,itot),wanall(itot) 1607 enddo ! itot 1608 iatomcor=0 1609 do itot=1,dimoverlap 1610 psichivect(ib,itot)=wanall(itot) 1611 enddo 1612 ! do iatom=1,natom 1613 ! if(paw_dmft%lpawu(iatom).ne.-1) then 1614 ! ndim=2*paw_dmft%lpawu(iatom)+1 1615 ! iatomcor=iatomcor+1 1616 ! do im=1,ndim 1617 ! do ispinor=1,nspinor 1618 ! paw_dmft%psichi(isppol,jkpt,ib,ispinor,iatom,im)=wanall(iatomcor,isppol,ispinor,im) 1619 ! end do ! ispinor 1620 ! end do ! im 1621 ! endif 1622 ! enddo ! iatom 1623 enddo ! ib 1624 1625 1626 ! calculation of overlap (check) 1627 largeoverlap=czero 1628 do ib=1,mbandc 1629 do itot=1,dimoverlap 1630 do itot1=1,dimoverlap 1631 largeoverlap(itot,itot1)=largeoverlap(itot,itot1)+ & 1632 & psichivect(ib,itot)*conjg(psichivect(ib,itot1)) 1633 enddo ! itot1 1634 enddo ! itot 1635 enddo ! ib 1636 1637 write(std_out,*)"jkpt=",jkpt 1638 do itot=1,dimoverlap 1639 write(std_out,'(100f7.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap) 1640 enddo 1641 1642 1643 ! psichivect -> psichi 1644 do ib=1,mbandc 1645 itot=0 1646 do iatom=1,natom 1647 if(paw_dmft%lpawu(iatom).ne.-1) then 1648 ndim=2*paw_dmft%lpawu(iatom)+1 1649 iatomcor=iatomcor+1 1650 do im=1,ndim 1651 do ispinor=1,nspinor 1652 itot=itot+1 1653 paw_dmft%psichi(isppol,jkpt,ib,ispinor,iatom,im)=psichivect(ib,itot) 1654 end do ! ispinor 1655 end do ! im 1656 endif 1657 enddo ! iatom 1658 enddo ! ib 1659 1660 1661 ! End big loop over isppol 1662 enddo !ispppol 1663 1664 ABI_FREE(psichivect) 1665 ABI_FREE(sqrtmatinv) 1666 ABI_FREE(wanall) 1667 ABI_FREE(largeoverlap) 1668 1669 endif 1670 1671 end subroutine normalizepsichi 1672 1673 end subroutine psichi_renormalization