TABLE OF CONTENTS
ABINIT/dmft_solve [ Functions ]
NAME
dmft_solve
FUNCTION
Solve the DMFT loop from PAW data.
INPUTS
cryst_struc <type(crystal_t)>=crystal structure data istep = step of iteration for DFT. dft_occup <type(oper_type)> = occupations in the correlated orbitals in DFT paw_dmft <type(paw_dmft_type)> = data for self-consistent DFT+DMFT calculations. pawang <type(pawang)>=paw angular mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data pawprtvol = option for printing
OUTPUT
paw_dmft <type(paw_dmft_type)> = data for self-consistent DFT+DMFT calculations.
NOTES
SOURCE
93 subroutine dmft_solve(cryst_struc,istep,dft_occup,paw_dmft,pawang,pawtab,pawprtvol) 94 95 !Arguments ------------------------------------ 96 !scalars 97 integer, intent(in) :: istep 98 integer, intent(in) :: pawprtvol 99 !type(MPI_type), intent(inout) :: mpi_enreg 100 type(pawang_type), intent(in) :: pawang 101 type(crystal_t),intent(in) :: cryst_struc 102 type(pawtab_type),intent(in) :: pawtab(cryst_struc%ntypat) 103 type(oper_type),intent(in) :: dft_occup 104 type(paw_dmft_type), intent(inout) :: paw_dmft 105 106 !Local variables ------------------------------ 107 !array 108 real(dp) :: tsec(2) 109 !scalars 110 integer :: check,idmftloop,istep_iter,spaceComm,my_rank,opt_renorm 111 integer :: itypat,natomcor,iatom 112 logical :: etot_var 113 character(len=200) :: char_enddmft 114 ! type 115 type(green_type) :: green 116 type(green_type) :: greendft 117 type(hu_type),allocatable :: hu(:) 118 type(green_type) :: weiss 119 type(self_type) :: self 120 type(self_type) :: self_new 121 !type(oper_type) :: self_minus_hdc_oper 122 type(energy_type) :: energies_dmft 123 type(energy_type) :: energies_tmp 124 character(len=500) :: message 125 character(len=5) :: thdyn 126 character(len=4) :: part2,part3 127 !************************************************************************ 128 129 DBG_ENTER('COLL') 130 my_rank = xmpi_comm_rank(paw_dmft%spacecomm) 131 132 check=paw_dmft%dmftcheck ! checks enabled 133 !paw_dmft%dmft_fermi_prec=tol5 134 paw_dmft%dmft_fermi_prec=paw_dmft%dmft_charge_prec*ten 135 !paw_dmft%dmft_charge_prec=20_dp ! total number of electron. 136 paw_dmft%dmft_prgn=1 137 paw_dmft%dmft_prgn=0 138 etot_var=.true. 139 thdyn="fcalc" 140 thdyn="ecalc" 141 if(thdyn=="ecalc") then ! valid 142 part2="both" 143 part3="none" 144 else if(thdyn=="fcalc") then ! not tested 145 part2="corr" 146 part3="band" 147 end if 148 149 if(check==1) then 150 write(message,'(2a)') ch10,' DMFT Checks are enabled ' 151 else 152 write(message,'(2a)') ch10,' DMFT Checks will not be performed' 153 end if 154 call wrtout(std_out,message,'COLL') 155 156 157 if(istep==0) then 158 message = ' istep should not be equal to zero' 159 ABI_BUG(message) 160 end if 161 spaceComm=paw_dmft%spacecomm 162 !if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt 163 !call xmpi_barrier(spaceComm) 164 165 call initialize_self(self,paw_dmft) 166 call init_energy(cryst_struc,energies_dmft) 167 168 !=========================================================================== 169 !== First construct DFT green function (Init, Compute, Integrate, Print) 170 !=========================================================================== 171 write(message,'(6a)') ch10,' =====================================', & 172 & ch10,' ===== DFT Green Function Calculation',& 173 & ch10,' =====================================' 174 call wrtout(std_out,message,'COLL') 175 call icip_green("DFT",cryst_struc,greendft,paw_dmft,pawang,3,self) 176 !call print_green('DFT_NOT_renormalized',greendft,1,paw_dmft,pawprtvol=1,opt_wt=1) 177 178 !== Compare greendft%occup and dft_occup: check that DFT green function is fine 179 !---------------------------------------------------------------------- 180 write(message,'(2a)') ch10,& 181 & ' == Check lda occ. mat. from green with respect to the direct calc ==' 182 call wrtout(std_out,message,'COLL') 183 call diff_oper("Occup from DFT green function",& 184 & "DFT occupations",greendft%occup,dft_occup,1,paw_dmft%dmft_tolfreq) 185 ! write(message,'(2a)') ch10,& 186 !& ' ***** => Warning : diff_oper is suppressed for test' 187 ! call wrtout(std_out,message,'COLL') 188 write(message,'(2a)') ch10,& 189 & ' ***** => Calculation of Green function is thus correct without self ****' 190 call wrtout(std_out,message,'COLL') 191 call destroy_green(greendft) 192 193 !== Orthonormalize psichi 194 !---------------------------------------------------------------------- 195 call timab(621,1,tsec) 196 natomcor=0 197 do iatom=1,paw_dmft%natom 198 if(paw_dmft%lpawu(iatom).ne.-1) then 199 natomcor=natomcor+1 200 end if 201 end do 202 opt_renorm=paw_dmft%dmft_wanorthnorm 203 !if(natomcor>1) opt_renorm=2 ! if number of atoms is larger than one, one must use a new orthogonalization scheme. 204 if(paw_dmft%nspinor==2.and.(paw_dmft%dmft_solv==8.or.paw_dmft%dmft_solv==9)) opt_renorm=2 ! necessary to use hybri_limit in qmc_prep_ctqmc 205 ! ought to be generalized in the future 206 if(paw_dmft%dmft_solv/=-1) then 207 call psichi_renormalization(cryst_struc,paw_dmft,pawang,opt=opt_renorm) 208 209 210 ! check that loc_oper(upfold_oper)=I 211 ! call init_oper(paw_dmft,self_minus_hdc_oper) 212 ! call identity_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom) 213 ! call print_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom,1) 214 ! call upfold_oper(self_minus_hdc_oper,paw_dmft,1) 215 ! call print_oper(self_minus_hdc_oper,9,paw_dmft,3) 216 ! call loc_oper(self_minus_hdc_oper,paw_dmft,1) 217 ! call sym_matlu(cryst_struc,self_minus_hdc_oper%matlu,pawang) 218 ! call print_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom,1) 219 ! call destroy_oper(self_minus_hdc_oper) 220 221 222 ! =========================================================================== 223 ! == re-construct DFT green function with new psichis 224 ! =========================================================================== 225 write(message,'(6a)')& 226 & ch10,' ==============================================================='& 227 & ,ch10,' ===== DFT Green Function Calculation with renormalized psichi'& 228 & ,ch10,' ==============================================================' 229 end if 230 call timab(621,2,tsec) 231 232 call wrtout(std_out,message,'COLL') 233 call icip_green("DFT_renormalized",cryst_struc,greendft,& 234 & paw_dmft,pawang,pawprtvol,self) 235 !call print_green('DFT_renormalized',greendft,1,paw_dmft,pawprtvol=1,opt_wt=1) 236 237 !Need to store idmftloop and set it to zero to avoid useless print_energy in ab_out 238 idmftloop=paw_dmft%idmftloop 239 paw_dmft%idmftloop=0 240 call compute_energy(cryst_struc,energies_dmft,greendft,paw_dmft,pawprtvol,pawtab,self,occ_type=" lda",part='both') 241 paw_dmft%idmftloop=idmftloop 242 243 if(paw_dmft%dmft_prgn==1) then 244 if(paw_dmft%lpsichiortho==1) then 245 call local_ks_green(greendft,paw_dmft,prtopt=1) 246 end if 247 end if 248 call destroy_self(self) 249 !call printocc_green(greendft,9,paw_dmft,3,chtype="DFT GREEN PSICHI") 250 251 write(message,'(6a)')& 252 & ch10,' ==============================================================='& 253 & ,ch10,' ===== Define Interaction and self-energy' & 254 & ,ch10,' ==============================================================' 255 call wrtout(std_out,message,'COLL') 256 !== define Interaction from input upawu and jpawu 257 !---------------------------------------------------------------------- 258 ABI_MALLOC(hu,(cryst_struc%ntypat)) 259 call init_hu(cryst_struc,pawtab,hu,paw_dmft%dmftqmc_t2g,paw_dmft%dmftqmc_x2my2d) 260 call initialize_self(self,paw_dmft) 261 262 ! Set Hu in density representation for calculation of entropy if needed... 263 do itypat=1,cryst_struc%ntypat 264 if ( hu(itypat)%lpawu /= -1 ) then 265 call data4entropyDMFT_setHu(paw_dmft%forentropyDMFT,itypat,hu(itypat)%udens(:,:)) 266 end if 267 end do 268 269 !== define self from scratch or file and double counting 270 !---------------------------------------------------------------------- 271 !- Self allocated 272 call dc_self(greendft%charge_matlu,cryst_struc,hu,self,paw_dmft%dmft_dc,pawprtvol) 273 274 !- Read self or do self=hdc 275 if(paw_dmft%dmft_solv==4) then 276 ! write(std_out,*) "shift before rw_self",self%qmc_shift(1) 277 call make_qmcshift_self(cryst_struc,hu,self) 278 end if 279 call timab(627,1,tsec) 280 call rw_self(self,paw_dmft,prtopt=2,opt_rw=1,istep_iter=1000*istep) 281 call timab(627,2,tsec) 282 283 !== If QMC is used, and self energy is read for file, then 284 !== one does NOT shifts the self-energy because it was already shifted when writed, 285 !== and thus then weiss will be shifted 286 !---------------------------------------------------------------------- 287 !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf==1) & 288 !& call make_qmcshift_self(cryst_struc,hu,self) 289 !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf/=1) & 290 !& call make_qmcshift_self(cryst_struc,hu,self,apply=.true.) 291 292 call destroy_green(greendft) ! destroy DFT green function 293 call print_self(self,"print_dc",paw_dmft,prtopt=2) 294 295 296 297 !=========================================================================== 298 !== construct green function with the self-energy. 299 !=========================================================================== 300 write(message,'(6a)') & 301 & ch10,' =================================================================', & 302 & ch10,' ===== Green Function Calculation with input self-energy ========', & 303 & ch10,' =================================================================' 304 call wrtout(std_out,message,'COLL') 305 call icip_green("Green_inputself",cryst_struc,green,& 306 & paw_dmft,pawang,pawprtvol,self,opt_self=1) 307 !call print_green('beforefermi_green',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 308 ! call abi_abort('COLL') 309 310 !== Find fermi level 311 !--------------------------------------------------------------------- 312 !write(message,'(2a,i3,13x,a)') ch10,' === Compute green function from self-energy' 313 call fermi_green(cryst_struc,green,paw_dmft,pawang,self) 314 !== define weiss field only for the local quantities (opt_oper=2) 315 !---------------------------------------------------------------------- 316 ! write(std_out,*) "nkpt befreo init_greenweiss",ifreq,paw_dmft%nkpt 317 call init_green(weiss,paw_dmft,opt_oper_ksloc=2) 318 ! do ifreq=1,weiss%nw 319 ! write(std_out,*) "nkpt from weiss1",ifreq,weiss%oper(ifreq)%nkpt 320 ! enddo 321 322 !== Check fourier transforms 323 !---------------------------------------------------------------------- 324 if(check==1) then 325 call check_fourier_green(cryst_struc,green,paw_dmft,pawang) 326 end if 327 328 !== If QMC is used, and self energy is not read for file, then 329 !== one shifts the self-energy, and thus then weiss will be shifted 330 !== after dyson, in a coherent way concerning qmc_shift and qmc_xmu. 331 !---------------------------------------------------------------------- 332 !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf/=1) & 333 !& call make_qmcshift_self(cryst_struc,hu,self,apply=.true.) 334 !if(paw_dmft%dmft_solv==4) write(std_out,*) "shift after make_qmcshift_self",self%qmc_shift(1) 335 336 write(message,'(6a)') & 337 & ch10,' ======================================================'& 338 & ,ch10,' ===== DMFT Loop starts here ========'& 339 & ,ch10,' ======================================================' 340 call wrtout(std_out,message,'COLL') 341 342 !======================================================================= 343 !=== dmft loop ======================================================= 344 do idmftloop=1, paw_dmft%dmft_iter 345 !paw_dmft%idmftloop=idmftloop 346 paw_dmft%idmftloop=paw_dmft%idmftloop+1 347 ! ======================================================================= 348 istep_iter=1000*istep+idmftloop 349 350 write(message,'(2a,i3,13x,a)') ch10,& 351 & ' ===== DMFT Loop : ITER number',paw_dmft%idmftloop,'========' 352 call wrtout(std_out,message,'COLL') 353 354 ! == Dyson Equation G,self -> weiss(w) 355 ! --------------------------------------------------------------------- 356 call dyson(green,paw_dmft,self,weiss,opt_weissself=1) 357 ! call print_green('afterDyson',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 358 ! call abi_abort('COLL') 359 360 ! == Printout local "occupations" from weiss field (useless) 361 if(abs(pawprtvol)>3) then 362 call integrate_green(cryst_struc,weiss,paw_dmft,& 363 & pawang,prtopt=2,opt_ksloc=2) 364 call printocc_green(weiss,5,paw_dmft,3,opt_weissgreen=1) 365 end if 366 367 ! === Prepare data, solve Impurity problem: weiss(w) -> G(w) 368 ! --------------------------------------------------------------------- 369 call initialize_self(self_new,paw_dmft) 370 371 call impurity_solve(cryst_struc,green,hu,& 372 & paw_dmft,pawang,pawtab,self,self_new,weiss,pawprtvol) ! weiss-> green, or self if dmft_solv=1 373 ! if(paw_dmft%dmft_solv==4) write(std_out,*) "shift after impurity",self%qmc_shift(1) 374 375 ! == Compute double counting from charge from green_solver 376 ! --------------------------------------------------------------------- 377 if (green%has_charge_matlu_solver/=2) green%charge_matlu_solver=green%charge_matlu 378 call dc_self(green%charge_matlu_solver,cryst_struc,hu,self_new,paw_dmft%dmft_dc,pawprtvol) 379 380 ! == Solve dyson equation. G_imp(w), weiss_imp(w) -> Self_imp(w) 381 ! --------------------------------------------------------------------- 382 ! if dmft_solv==1, self is computed previously 383 if(abs(paw_dmft%dmft_solv)/=1) then 384 call dyson(green,paw_dmft,self_new,weiss,opt_weissself=2) 385 end if 386 ! do ifreq=1,green%nw 387 ! call sym_matlu(cryst_struc,self%oper(ifreq)%matlu,pawang) 388 ! enddo 389 390 ! == Possibility if imposing self (opt_rw==3) 391 ! --------------------------------------------------------------------- 392 call timab(627,1,tsec) 393 call rw_self(self_new,paw_dmft,prtopt=2,opt_rw=3,istep_iter=istep_iter) 394 call timab(627,2,tsec) 395 396 ! print dc just computed before and self computed in before in dyson or 397 ! impurity_solve 398 if(abs(pawprtvol)>=3) then 399 write(message,'(2a)') ch10," == New self" 400 call wrtout(std_out,message,'COLL') 401 call print_self(self_new,"print_dc",paw_dmft,2) 402 write(message,'(2a)') ch10," == Old self" 403 call wrtout(std_out,message,'COLL') 404 call print_self(self,"print_dc",paw_dmft,2) 405 end if 406 407 ! if(paw_dmft%dmft_solv==4) write(std_out,*) "shift before computeenergy ",self%qmc_shift(1) 408 ! == Compute Energy with NEW self-energy and edc from green_solver, 409 ! new local green function and old occupations for eband 410 ! fermi level not optimized for this self_energy. 411 ! --------------------------------------------------------------------- 412 ! green= local green function and local charge comes directly from solver 413 ! green= ks green function and occupations comes from old_self 414 call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,& 415 & pawtab,self_new,occ_type="nlda",part=part2) 416 417 ! == Mix new and old self_energies and double countings 418 ! --------------------------------------------------------------------- 419 call new_self(self,self_new,paw_dmft,1) ! self,self_new => self 420 write(message,'(2a)') ch10," == After mixing," 421 !print *, " my_rank newself", my_rank,self%oper(1)%matlu(1)%mat(1,1,1,1,1) 422 call wrtout(std_out,message,'COLL') 423 call print_self(self,"print_dc",paw_dmft,2) ! print self and DC 424 call destroy_self(self_new) 425 426 ! == Compute green function self -> G(k) 427 ! --------------------------------------------------------------------- 428 call compute_green(cryst_struc,green,paw_dmft,pawang,1,self,opt_self=1,opt_nonxsum=1) 429 430 call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=3,opt_ksloc=3,opt_diff=1) !,opt_nonxsum=1) 431 432 call printocc_green(green,5,paw_dmft,3,chtype="DMFT") 433 ! call printocc_green(green,9,paw_dmft,3,chtype="DMFT FULL") 434 if(paw_dmft%lpsichiortho==1.and.paw_dmft%dmft_prgn==1) then 435 call local_ks_green(green,paw_dmft,prtopt=1) 436 end if 437 438 ! == Find fermi level 439 ! --------------------------------------------------------------------- 440 call fermi_green(cryst_struc,green,paw_dmft,pawang,self) 441 ! call abi_abort('COLL') 442 443 ! == Compute Energy with Mixed self-energy and green function recomputed with new self 444 ! --------------------------------------------------------------------- 445 ! green= lattice green function computed from self for a given chemical potential mu (self_mixed,mu) 446 ! green= local green function is computed from lattice green function(self_mixed,mu) 447 ! green= occupations are computed with lattice green function(self_mixed,mu) 448 call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,occ_type="nlda",part=part3) 449 450 ! == Save self on disk 451 ! --------------------------------------------------------------------- 452 call timab(627,1,tsec) 453 call rw_self(self,paw_dmft,prtopt=2,opt_rw=2,pawang=pawang,cryst_struc=cryst_struc) 454 call timab(627,2,tsec) 455 456 ! == Test convergency 457 ! --------------------------------------------------------------------- 458 char_enddmft="DMFT (end of DMFT loop)" 459 if(green%ifermie_cv==1.and.self%iself_cv==1.and.green%ichargeloc_cv==1.and.paw_dmft%idmftloop>1) then 460 write(message,'(a,8x,a)') ch10,"DMFT Loop is converged !" 461 call wrtout(std_out,message,'COLL') 462 char_enddmft="converged DMFT" 463 exit 464 end if 465 ! ======================================================================= 466 ! === end dmft loop ==================================================== 467 end do 468 !========================================================================= 469 470 !== Save self on disk 471 !------------------------------------------------------------------------- 472 call timab(627,1,tsec) 473 call rw_self(self,paw_dmft,prtopt=2,opt_rw=2) 474 call timab(627,2,tsec) 475 476 !paw_dmft%idmftloop=0 477 478 write(message,'(2a,13x,a)') ch10,' ===== DMFT Loop : END ',& 479 & '========' 480 call wrtout(std_out,message,'COLL') 481 482 if(paw_dmft%dmft_entropy/=0) then 483 ! compute Edc for U=1 and J=U/J 484 call init_energy(cryst_struc,energies_tmp) 485 !call compute_dftu_energy(cryst_struc,energies_tmp,green,paw_dmft,pawtab) 486 call compute_dftu_energy(cryst_struc,energies_tmp,green,paw_dmft,pawtab,paw_dmft%forentropyDMFT%J_over_U) 487 call data4entropyDMFT_setDc(paw_dmft%forentropyDMFT,energies_tmp%e_dc(:)) 488 call destroy_energy(energies_tmp,paw_dmft) 489 endif 490 491 !== Compute final values for green functions, occupations, and spectral function 492 !-------------------------------------------------------------------------------- 493 !Do not compute here, because, one want a energy computed after the 494 !solver (for Hubbard I and DFT+U). 495 call compute_green(cryst_struc,green,paw_dmft,pawang,1,self,opt_self=1,opt_nonxsum=1) 496 call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=2,opt_ksloc=3) !,opt_nonxsum=1) 497 !call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,opt=0) 498 idmftloop=paw_dmft%idmftloop 499 paw_dmft%idmftloop=0 500 call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,occ_type="nlda",part="band") 501 paw_dmft%idmftloop=idmftloop 502 503 !write(message,'(2a,13x,a)') ch10,' ===== DMFT Loop is finished' 504 !call wrtout(ab_out,message,'COLL') 505 !write(std_out,*) "PRINTOCC INITIAL" 506 call printocc_green(green,9,paw_dmft,3,chtype=char_enddmft) 507 !write(std_out,*) "KS=czero" 508 !green%occup%ks=czero 509 !write(std_out,*) "PRINTOCC AFTER KS=0" 510 !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT") 511 !write(std_out,*) "UPFOLD_OPER" 512 !call upfold_oper(green%occup,paw_dmft,1) 513 !write(std_out,*) "PRINTOCC AFTER UPFOLD_OPER" 514 !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT") 515 !write(std_out,*) "MATLU=czero" 516 !green%occup%matlu(1)%mat=czero 517 !green%occup%ks(:,:,:,:)=cmplx(real(green%occup%ks(:,:,:,:))) 518 !write(std_out,*) "PRINTOCC AFTER MATLU=0 AND IMAG KS=0" 519 !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT") 520 !write(std_out,*) "LOC_OPER" 521 !call loc_oper(green%occup,paw_dmft,1) 522 !write(std_out,*) "PRINTOCC AFTER LOC_OPER" 523 !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT") 524 !call flush_unit(std_out) 525 !call abi_abort('COLL') 526 if(paw_dmft%dmft_solv<=2.and.paw_dmft%prtdos>=1) then 527 call spectral_function(cryst_struc,green,hu,paw_dmft,pawang,pawtab,self,pawprtvol) 528 end if 529 call destroy_green(weiss) 530 call destroy_green(green) 531 !todo_ab rotate back density matrix into unnormalized basis just for 532 !printout 533 call destroy_hu(hu,cryst_struc%ntypat,paw_dmft%dmftqmc_t2g,paw_dmft%dmftqmc_x2my2d) 534 call destroy_self(self) 535 call destroy_energy(energies_dmft,paw_dmft) 536 537 write(message,'(2a,13x,a)') ch10,' ===== DMFT : END ',& 538 & '========' 539 call wrtout(std_out,message,'COLL') 540 541 ABI_FREE(hu) 542 543 DBG_EXIT("COLL") 544 545 end subroutine dmft_solve
ABINIT/dyson [ Functions ]
NAME
dyson
FUNCTION
Use the Dyson Equation to compute self-energy from green function
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset istep = step of iteration for DFT. dft_occup mpi_enreg=information about MPI parallelization paw_dmft = data for self-consistent DFT+DMFT calculations. opt_weissgreen= 1: compute weiss from green and self = 2: compute green from weiss and self = 4: compute green from weiss and self without inversion of weiss (weiss is previously inverted)
OUTPUT
edmft = energy in DMFT. paw_dmft = data for self-consistent DFT+DMFT calculations.
NOTES
SOURCE
828 subroutine dyson(green,paw_dmft,self,weiss,opt_weissself) 829 830 !Arguments ------------------------------------ 831 !scalars 832 type(green_type),intent(inout) :: green 833 type(paw_dmft_type), intent(inout) :: paw_dmft 834 type(self_type),intent(inout) :: self 835 type(green_type),intent(inout) :: weiss 836 integer,intent(in) :: opt_weissself 837 ! type(paw_dmft_type), intent(inout) :: paw_dmft 838 839 !Local variables ------------------------------ 840 !arrays 841 real(dp) :: tsec(2) 842 !scalars 843 integer :: ifreq,natom,nspinor,nsppol,weissinv 844 type(green_type) :: greeninv 845 character(len=500) :: message 846 ! type 847 ! type(matlu_type), pointer :: matlutemp,matlu1,matlu2 848 !************************************************************************ 849 850 call timab(623,1,tsec) 851 DBG_ENTER("COLL") 852 natom=green%oper(1)%natom 853 nsppol=green%oper(1)%nsppol 854 nspinor=green%oper(1)%nspinor 855 weissinv=1 856 if(opt_weissself==1) then 857 write(message,'(2a,i3,13x,a)') ch10,' === Use Dyson Equation => weiss ' 858 call wrtout(std_out,message,'COLL') 859 else if(opt_weissself==2) then 860 write(message,'(2a,i3,13x,a)') ch10,' === Use Dyson Equation => self' 861 call wrtout(std_out,message,'COLL') 862 end if 863 864 !call xmpi_barrier(spaceComm) 865 if(paw_dmft%dmft_solv==2) weissinv=0 866 call init_green(greeninv,paw_dmft,opt_oper_ksloc=2,wtype=green%w_type) 867 call copy_green(green,greeninv,opt_tw=2) 868 869 do ifreq=1,green%nw 870 if(opt_weissself==1) then 871 call inverse_oper(greeninv%oper(ifreq),option=1,prtopt=1) 872 ! warning green is now inversed 873 call add_matlu(greeninv%oper(ifreq)%matlu,self%oper(ifreq)%matlu,& 874 & weiss%oper(ifreq)%matlu,natom,sign_matlu2=1) 875 call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1) 876 else if(opt_weissself==2) then 877 878 ! write(59,*) paw_dmft%omega_lo(ifreq), real(green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)) 879 ! write(61,*) paw_dmft%omega_lo(ifreq), real(greeninv%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(greeninv%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)) 880 ! write(60,*) paw_dmft%omega_lo(ifreq), real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)) 881 ! call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1) 882 ! write(62,*) paw_dmft%omega_lo(ifreq), real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)) 883 ! call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1) 884 ! write(63,*) paw_dmft%omega_lo(ifreq), real(self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)) 885 886 ! write(std_out,*) "-----------------------IFREQ",ifreq 887 ! call print_matlu(greeninv%oper(ifreq)%matlu,paw_dmft%natom,1,opt_diag=-1) 888 call inverse_oper(greeninv%oper(ifreq),option=1,prtopt=1) 889 ! call print_matlu(greeninv%oper(ifreq)%matlu,paw_dmft%natom,1,opt_diag=-1) 890 ! If paw_dmft%dmft_solv==2, then inverse of weiss function is 891 ! computed in m_hubbard_one.F90 892 if(weissinv/=0) then 893 call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1) 894 end if 895 896 ! write(std_out,*) weiss%oper(1)%matlu(ifreq)%mat(1,1,1,1,1),"-",greeninv%oper(ifreq) 897 call add_matlu(weiss%oper(ifreq)%matlu,greeninv%oper(ifreq)%matlu,& 898 & self%oper(ifreq)%matlu,natom,sign_matlu2=-1) 899 900 ! write(64,*) paw_dmft%omega_lo(ifreq), real(greeninv%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(greeninv%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)) 901 ! write(65,*) paw_dmft%omega_lo(ifreq), real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)) 902 ! write(66,*) paw_dmft%omega_lo(ifreq), real(self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),aimag(self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)) 903 else 904 message = " BUG in dyson.F90" 905 ABI_BUG(message) 906 end if 907 end do 908 909 call destroy_green(greeninv) 910 911 912 call timab(623,2,tsec) 913 DBG_EXIT("COLL") 914 915 end subroutine dyson
ABINIT/impurity_solve [ Functions ]
NAME
impurity_solve
FUNCTION
Solve the Impurity problem
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
570 subroutine impurity_solve(cryst_struc,green,hu,paw_dmft,& 571 & pawang,pawtab,self_old,self_new,weiss,pawprtvol) 572 573 !Arguments ------------------------------------ 574 !scalars 575 ! type(pawang_type), intent(in) :: pawang 576 type(crystal_t),intent(in) :: cryst_struc 577 type(green_type), intent(inout) :: weiss 578 type(green_type), intent(inout) :: green !vz_i 579 type(hu_type),intent(inout) :: hu(cryst_struc%ntypat) 580 !type(MPI_type), intent(in) :: mpi_enreg 581 type(pawang_type), intent(in) :: pawang 582 type(pawtab_type),intent(in) :: pawtab(cryst_struc%ntypat) 583 type(paw_dmft_type), intent(inout) :: paw_dmft 584 type(self_type), intent(inout) :: self_new 585 type(self_type), intent(inout) :: self_old 586 integer, intent(in) :: pawprtvol 587 588 !Local variables ------------------------------ 589 real(dp) :: tsec(2) 590 character(len=500) :: message 591 complex(dpc) :: xx 592 integer :: ifreq 593 ! integer iatom,il,i_nd,isppol,lpawu,im,Nd,nrat,nsweeptot 594 ! real(dp) :: acc,kx 595 ! real(dp), allocatable :: correl(:,:),g0(:,:),gtmp(:,:) 596 !scalars 597 !************************************************************************ 598 !character(len=500) :: message 599 600 call timab(622,1,tsec) 601 !======================================================================= 602 !== Prepare data for Hirsch Fye QMC 603 !== NB: for CTQMC, Fourier Transformation are done inside the CTQMC code 604 !======================================================================= 605 if(abs(paw_dmft%dmft_solv)==4) then 606 ! == Initialize weiss and green functions for fourier transformation 607 ! ------------------------------------------------------------------- 608 write(message,'(2a,i3,13x,a)') ch10,' === Initialize Weiss field G_0(tau)' 609 call wrtout(std_out,message,'COLL') 610 call init_green_tau(weiss,paw_dmft) 611 call init_green_tau(green,paw_dmft) 612 ! in init_solver 613 614 ! == Print weiss function G_0(tau=0-) before computation (really useless check) 615 ! ------------------------------------------------------------------------------ 616 if(abs(pawprtvol)>3) then 617 write(message,'(2a,i3,13x,a)') ch10,' === Check G_0(tau=0-) first' 618 call wrtout(std_out,message,'COLL') 619 call printocc_green(weiss,6,paw_dmft,3) 620 end if 621 622 ! == Fourier transform of weiss Field 623 ! ------------------------------------ 624 ! for fourier of KS green functions 625 ! call fourier_green(cryst_struc,weiss,mpi_enreg,paw_dmft,pawang,pawtab,1) 626 write(message,'(2a,i3,13x,a)') ch10,' === Inverse Fourier Transform w->t of Weiss Field' 627 call wrtout(std_out,message,'COLL') 628 call fourier_green(cryst_struc,weiss,paw_dmft,pawang,opt_ksloc=2,opt_tw=-1) 629 630 ! == Print weiss function G2_0(tau=0-) 631 ! -------------------------------------- 632 call printocc_green(weiss,6,paw_dmft,3,opt_weissgreen=1) 633 634 ! for fourier of KS green functions 635 ! call fourier_green(cryst_struc,weiss,mpi_enreg,paw_dmft,pawang,pawtab,1) 636 ! == Print G_0(tau) in files 637 ! --------------------------- 638 if(paw_dmft%dmft_prgn==1) then 639 call print_green('weiss',weiss,1,paw_dmft,pawprtvol=1,opt_wt=2) 640 end if 641 642 else if(abs(paw_dmft%dmft_solv)>=5) then 643 ! == Initialize green functions for imaginary times 644 ! ------------------------------------------------------------------- 645 write(message,'(2a,i3,13x,a)') ch10,' === Initialize Green function G(tau)' 646 call wrtout(std_out,message,'COLL') 647 call init_green_tau(green,paw_dmft) 648 649 end if 650 !======================================================================= 651 !== End preparation of QMC 652 !======================================================================= 653 654 !======================================================================= 655 !== Solve impurity model ============================================= 656 !======================================================================= 657 write(message,'(2a,i3,13x,a)') ch10,' === Solve impurity model' 658 call wrtout(std_out,message,'COLL') 659 if(abs(paw_dmft%dmft_solv)==1) then 660 661 ! == DFT+U for test -> self 662 ! ------------------- 663 call dftu_self(cryst_struc,green,paw_dmft,& 664 & pawtab,self_new,opt_dftu=1,prtopt=pawprtvol) 665 else if(abs(paw_dmft%dmft_solv)==2) then 666 667 ! == Hubbard One -> green 668 ! ------------------- 669 call hubbard_one(cryst_struc,green,hu,paw_dmft,& 670 & pawang,pawprtvol,self_old%hdc,weiss) 671 672 else if(abs(paw_dmft%dmft_solv)==4) then 673 674 ! == QMC 675 ! ------------------- 676 call copy_green(weiss,green,opt_tw=1) 677 ! call qmc_prep 678 message = ' === QMC not yet distributed ' 679 ABI_ERROR(message) 680 ! call qmc_prep(cryst_struc,green,hu,mpi_enreg,paw_dmft,pawang& 681 !& ,pawprtvol,self_old%qmc_xmu,weiss) 682 683 else if(abs(paw_dmft%dmft_solv)>=5) then 684 685 ! == Nothing 686 ! ------------------- 687 ! call copy_green(weiss,green,opt_tw=1) 688 ! call copy_green(weiss,green,opt_tw=2) 689 690 call qmc_prep_ctqmc(cryst_struc,green,self_old,hu,paw_dmft,pawang,pawprtvol,weiss) 691 692 693 else if(abs(paw_dmft%dmft_solv)==0) then 694 695 ! == Nothing 696 ! ------------------- 697 ! weiss%occup%has_operks=0 -> only local part is duplicated 698 call copy_green(weiss,green,opt_tw=2) 699 end if 700 !call print_green("invWeiss",cryst_struc,weiss,3,paw_dmft,pawtab,2) 701 702 !======================================================================= 703 !== Treat data from HF QMC 704 !======================================================================= 705 if(abs(paw_dmft%dmft_solv)>=4) then 706 ! propagate qmc_shift (useful for compute_energy) 707 if(abs(paw_dmft%dmft_solv)==4) then 708 self_new%qmc_shift(:)=self_old%qmc_shift(:) 709 self_new%qmc_xmu(:)=self_old%qmc_xmu(:) 710 end if 711 712 ! == Print local occupations from G(tau) 713 ! --------------------------------------- 714 715 ! == Fourier back transform of green function G(tau)->G(iw_n) and 716 ! == compute occupations from g(tau) 717 ! ------------------------------------------------------------------- 718 if(abs(paw_dmft%dmft_solv)==4) then 719 write(message,'(2a,i3,13x,a)') ch10,' === Direct Fourier Transform t->w of Green Function' 720 call wrtout(std_out,message,'COLL') 721 call fourier_green(cryst_struc,green,paw_dmft,& 722 & pawang,opt_ksloc=2,opt_tw=1) 723 do ifreq=1,green%nw 724 xx= green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) 725 write(112,*) paw_dmft%omega_lo(ifreq),real(one/xx),aimag(one/xx) 726 write(113,*) paw_dmft%omega_lo(ifreq),real(xx),aimag(xx) 727 end do 728 call flush_unit(112) 729 call flush_unit(113) 730 ! if(paw_dmft%dmft_solv==5) stop 731 if(pawprtvol>=3) then 732 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 733 & " == Print green function for small freq after fourier " ! debug 734 call wrtout(std_out,message,'COLL') ! debug 735 call print_matlu(green%oper(1)%matlu,paw_dmft%natom,1) ! debug 736 end if 737 738 write(message,'(2a,i3,13x,a)') ch10,' INVERSE FOURIER OF G0 SUPPRESSED' 739 call wrtout(std_out,message,'COLL') 740 end if 741 if(abs(paw_dmft%dmft_solv)==888) then 742 ! == Back fourier transform of G_0(tau) for compensation (try or comment or improve FT). 743 ! ------------------------------------------------------------------- 744 write(message,'(2a,i3,13x,a)') ch10,' === Direct Fourier transform t->w of Weiss' 745 call wrtout(std_out,message,'COLL') 746 call fourier_green(cryst_struc,weiss,paw_dmft,& 747 & pawang,opt_ksloc=2,opt_tw=1) 748 749 if(pawprtvol>=3) then 750 write(message,'(a,2x,a,f13.5)') ch10,& ! debug 751 & " == Print weiss function for small freq after fourier " ! debug 752 call wrtout(std_out,message,'COLL') ! debug 753 call print_matlu(weiss%oper(1)%matlu,paw_dmft%natom,1) ! debug 754 end if 755 call destroy_green_tau(weiss) 756 end if 757 758 ! == Destroy tau part of green 759 ! ------------------------------------------------------------------- 760 call trace_oper(green%occup_tau,green%charge_ks,green%charge_matlu_solver,2) 761 green%has_charge_matlu_solver=2 762 call destroy_green_tau(green) 763 end if 764 !======================================================================= 765 !== End Treat data for QMC 766 !======================================================================= 767 768 !======================================================================= 769 !== Integrate green function and printout occupations 770 !======================================================================= 771 !For dmft_solv=-1,0,or 1 , the green function was not yet computed: it 772 !cannot be integrated 773 !======================================================================= 774 if(paw_dmft%dmft_solv>=2.and.green%w_type=="imag") then 775 ! == Integrate G(iw_n) 776 ! --------------------- 777 write(message,'(2a,i3,13x,a)') ch10,' === Integrate local part of green function' 778 call wrtout(std_out,message,'COLL') 779 call integrate_green(cryst_struc,green,paw_dmft,& 780 & pawang,prtopt=2,opt_ksloc=2,opt_after_solver=1) 781 782 ! == Print local occupations from integration of G(iw_n) 783 ! -------------------------------------------------------- 784 call printocc_green(green,5,paw_dmft,3) 785 786 ! == Print G_loc(w) 787 ! -------------------------------------------------------- 788 if(paw_dmft%dmft_prgn==1) then 789 call print_green('DMFT_IMPURITY',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 790 end if 791 end if 792 !stop 793 794 if(abs(pawprtvol)>0) then 795 end if 796 797 call timab(622,2,tsec) 798 end subroutine impurity_solve
ABINIT/m_dmft [ Modules ]
NAME
m_dmft
FUNCTION
COPYRIGHT
Copyright (C) 2006-2022 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
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 24 #include "abi_common.h" 25 26 MODULE m_dmft 27 28 use defs_basis 29 #ifdef HAVE_NETCDF 30 use netcdf 31 #endif 32 use m_xmpi 33 use m_errors 34 use m_abicore 35 use m_data4entropyDMFT 36 37 use m_time, only : timab 38 use m_pawang, only : pawang_type 39 use m_pawtab, only : pawtab_type 40 use m_paw_dmft, only: paw_dmft_type 41 use m_crystal, only : crystal_t 42 use m_green, only : green_type, destroy_green, icip_green,init_green,& 43 & print_green,printocc_green,& 44 & integrate_green,copy_green,compute_green,& 45 & check_fourier_green,local_ks_green,fermi_green, & 46 & fourier_green, init_green_tau,destroy_green_tau 47 use m_oper, only : oper_type,diff_oper,upfold_oper,loc_oper, trace_oper, inverse_oper !,upfold_oper,init_oper,destroy_oper,print_oper 48 use m_self, only : self_type,initialize_self,destroy_self,print_self,dc_self,rw_self,new_self,make_qmcshift_self 49 use m_hu, only : hu_type,init_hu,destroy_hu 50 use m_energy, only : energy_type,init_energy,destroy_energy,compute_energy,print_energy,compute_dftu_energy 51 use m_matlu, only : print_matlu,sym_matlu, matlu_type,init_matlu,destroy_matlu, add_matlu, copy_matlu 52 use m_datafordmft, only : psichi_renormalization 53 use m_io_tools, only : flush_unit 54 use m_hubbard_one, only : hubbard_one 55 use m_dftu_self, only : dftu_self 56 use m_forctqmc, only : qmc_prep_ctqmc 57 58 implicit none 59 60 private 61 62 public :: dmft_solve 63 public :: impurity_solve 64 public :: dyson 65 public :: spectral_function
m_dmft/spectral_function [ Functions ]
[ Top ] [ m_dmft ] [ Functions ]
NAME
spectral_function
FUNCTION
Print the spectral function computed from Green function in real frequency
INPUTS
cryst_struc <type(crystal_t)>=crystal structure data green <type(green_type)>= green function data hu <type(hu_type)>= datatype of type hu paw_dmft = data for self-consistent DFT+DMFT calculations. pawang <type(pawang)>=paw angular mesh and related data self <type(self_type)>= variables related to self-energy prtopt= option for printing
OUTPUT
paw_dmft = data for self-consistent DFT+DMFT calculations.
NOTES
SOURCE
941 subroutine spectral_function(cryst_struc,green,hu,paw_dmft,& 942 & pawang,pawtab,self_old,prtopt) 943 944 !Arguments ------------------------------------ 945 !scalars 946 ! type(pawang_type), intent(in) :: pawang 947 type(crystal_t),intent(in) :: cryst_struc 948 type(green_type), intent(in) :: green 949 type(hu_type),intent(inout) :: hu(cryst_struc%ntypat) 950 !type(MPI_type), intent(inout) :: mpi_enreg 951 type(pawang_type), intent(in) :: pawang 952 type(pawtab_type),intent(in) :: pawtab(cryst_struc%ntypat) 953 type(self_type), intent(inout) :: self_old 954 type(paw_dmft_type), intent(inout) :: paw_dmft 955 integer, intent(in) :: prtopt 956 957 !Local variables ------------------------------ 958 character(len=500) :: message 959 type(green_type) :: greenr 960 type(green_type) :: weissr 961 type(self_type) :: selfr 962 !scalars 963 !************************************************************************ 964 !character(len=500) :: message 965 966 ! opt_oper_ksloc=3 to be able to compute spectral function. 967 call init_green(greenr,paw_dmft,opt_oper_ksloc=3,wtype="real") 968 call init_green(weissr,paw_dmft,wtype="real") 969 call copy_matlu(green%occup%matlu,greenr%occup%matlu,paw_dmft%natom) 970 call initialize_self(selfr,paw_dmft,wtype="real") 971 !======================================================================= 972 !== Solve impurity model with green function for real frequency 973 !======================================================================= 974 write(message,'(2a,i3,13x,a)') ch10,' === Write Spectral function' 975 call wrtout(std_out,message,'COLL') 976 if(abs(paw_dmft%dmft_solv)==1) then 977 978 ! == DFT+U for test 979 ! ------------------- 980 call dftu_self(cryst_struc,greenr,paw_dmft,& 981 & pawtab,selfr,opt_dftu=1,prtopt=prtopt) 982 else if(abs(paw_dmft%dmft_solv)==2) then 983 984 ! == Hubbard One 985 ! ------------------- 986 call hubbard_one(cryst_struc,greenr,hu,paw_dmft,& 987 & pawang,prtopt,self_old%hdc,weissr) 988 989 else if(abs(paw_dmft%dmft_solv)==4) then 990 991 ! == Nothing 992 ! ------------------- 993 message = "spectral_function: This section of code is disabled!" 994 ABI_ERROR(message) 995 call copy_green(weissr,greenr,opt_tw=1) 996 997 else if(abs(paw_dmft%dmft_solv)>=5) then 998 999 ! == Nothing 1000 ! ------------------- 1001 ABI_ERROR("Stopping before copy_green") 1002 call copy_green(weissr,greenr,opt_tw=1) 1003 1004 else if(abs(paw_dmft%dmft_solv)==0) then 1005 1006 ! == Nothing 1007 ! ------------------- 1008 ! weiss%occup%has_operks=0 -> only local part is duplicated 1009 call copy_green(weissr,greenr,opt_tw=2) 1010 end if 1011 1012 !======================================================================= 1013 !== Integrate green function and printout occupations 1014 !For dmft_solv=-1,0,or 1 , the green function was not computed: it 1015 !cannot be integrated 1016 !======================================================================= 1017 call dc_self(green%charge_matlu_solver,cryst_struc,hu,selfr,paw_dmft%dmft_dc,prtopt) 1018 if(abs(paw_dmft%dmft_solv)/=1.and.paw_dmft%dmft_solv/=0) then 1019 call dyson(greenr,paw_dmft,selfr,weissr,opt_weissself=2) 1020 end if 1021 call compute_green(cryst_struc,greenr,paw_dmft,pawang,1,selfr,opt_self=1) 1022 call print_green("realw",greenr,4,paw_dmft,pawprtvol=3) 1023 call rw_self(selfr,paw_dmft,prtopt=2,opt_rw=2) 1024 1025 if(abs(prtopt)>0) then 1026 end if 1027 call destroy_self(selfr) 1028 call destroy_green(weissr) 1029 call destroy_green(greenr) 1030 1031 end subroutine spectral_function