TABLE OF CONTENTS
ABINIT/dmft_solve [ Functions ]
NAME
dmft_solve
FUNCTION
Solve the DMFT loop from PAW data.
COPYRIGHT
Copyright (C) 1999-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cryst_struc <type(crystal_t)>=crystal structure data istep = step of iteration for LDA. lda_occup <type(oper_type)> = occupations in the correlated orbitals in LDA paw_dmft <type(paw_dmft_type)> = data for self-consistent LDA+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 LDA+DMFT calculations.
NOTES
PARENTS
vtorho
CHILDREN
check_fourier_green,compute_energy,compute_green,compute_ldau_energy data4entropydmft_setdc,data4entropydmft_sethu,dc_self,destroy_energy destroy_green,destroy_hu,destroy_self,diff_oper,dyson,fermi_green icip_green,impurity_solve,init_energy,init_green,init_hu initialize_self,integrate_green,local_ks_green,make_qmcshift_self new_self,print_self,printocc_green,psichi_renormalization,rw_self spectral_function,timab,wrtout
SOURCE
44 #if defined HAVE_CONFIG_H 45 #include "config.h" 46 #endif 47 48 49 #include "abi_common.h" 50 51 subroutine dmft_solve(cryst_struc,istep,lda_occup,paw_dmft,pawang,pawtab,pawprtvol) 52 53 54 use defs_basis 55 use defs_abitypes 56 use m_xmpi 57 use m_errors 58 use m_profiling_abi 59 60 use m_pawang, only : pawang_type 61 use m_pawtab, only : pawtab_type 62 use m_paw_dmft, only: paw_dmft_type 63 use m_crystal, only : crystal_t 64 use m_green, only : green_type, destroy_green, icip_green,init_green,& 65 & print_green,printocc_green,& 66 & integrate_green,copy_green,compute_green,check_fourier_green 67 use m_oper, only : oper_type,diff_oper,upfold_oper,loc_oper 68 use m_self, only : self_type,initialize_self,destroy_self,print_self,dc_self,rw_self,new_self,make_qmcshift_self 69 use m_hu, only : hu_type,init_hu,destroy_hu 70 use m_energy, only : energy_type,init_energy,destroy_energy,compute_energy,print_energy,compute_ldau_energy 71 use m_matlu, only : print_matlu,sym_matlu 72 use m_data4entropyDMFT 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'dmft_solve' 78 use interfaces_14_hidewrite 79 use interfaces_18_timing 80 use interfaces_68_dmft, except_this_one => dmft_solve 81 !End of the abilint section 82 83 implicit none 84 85 !Arguments ------------------------------------ 86 !scalars 87 integer, intent(in) :: istep 88 integer, intent(in) :: pawprtvol 89 !type(MPI_type), intent(inout) :: mpi_enreg 90 type(pawang_type), intent(in) :: pawang 91 type(crystal_t),intent(in) :: cryst_struc 92 type(pawtab_type),intent(in) :: pawtab(cryst_struc%ntypat) 93 type(oper_type),intent(in) :: lda_occup 94 type(paw_dmft_type), intent(inout) :: paw_dmft 95 96 !Local variables ------------------------------ 97 !array 98 real(dp) :: tsec(2) 99 !scalars 100 integer :: check,idmftloop,istep_iter,spaceComm,my_rank,opt_renorm 101 integer :: itypat 102 logical :: etot_var 103 character(len=200) :: char_enddmft 104 ! type 105 type(green_type) :: green 106 type(green_type) :: greenlda 107 type(hu_type),allocatable :: hu(:) 108 type(green_type) :: weiss 109 type(self_type) :: self 110 type(self_type) :: self_new 111 type(energy_type) :: energies_dmft 112 type(energy_type) :: energies_tmp 113 character(len=500) :: message 114 character(len=5) :: thdyn 115 character(len=4) :: part2,part3 116 !************************************************************************ 117 118 DBG_ENTER('COLL') 119 my_rank = xmpi_comm_rank(paw_dmft%spacecomm) 120 121 check=paw_dmft%dmftcheck ! checks enabled 122 paw_dmft%dmft_fepr=tol5 123 paw_dmft%dmft_chpr=tol6 124 !paw_dmft%dmft_chpr=20_dp ! total number of electron. 125 paw_dmft%dmft_prgn=1 126 paw_dmft%dmft_prgn=0 127 etot_var=.true. 128 thdyn="fcalc" 129 thdyn="ecalc" 130 if(thdyn=="ecalc") then ! valid 131 part2="both" 132 part3="none" 133 else if(thdyn=="fcalc") then ! not tested 134 part2="corr" 135 part3="band" 136 end if 137 138 if(check==1) then 139 write(message,'(2a)') ch10,' DMFT Checks are enabled ' 140 else 141 write(message,'(2a)') ch10,' DMFT Checks will not be performed' 142 end if 143 call wrtout(std_out,message,'COLL') 144 145 146 if(istep==0) then 147 message = ' istep should not be equal to zero' 148 MSG_BUG(message) 149 end if 150 spaceComm=paw_dmft%spacecomm 151 !if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt 152 !call xmpi_barrier(spaceComm) 153 154 call initialize_self(self,paw_dmft) 155 call init_energy(cryst_struc,energies_dmft) 156 157 !=========================================================================== 158 !== First construct LDA green function (Init, Compute, Integrate, Print) 159 !=========================================================================== 160 write(message,'(6a)') ch10,' =====================================', & 161 & ch10,' ===== LDA Green Function Calculation',& 162 & ch10,' =====================================' 163 call wrtout(std_out,message,'COLL') 164 call icip_green("LDA",cryst_struc,greenlda,paw_dmft,pawang,3,self) 165 !call print_green('LDA_NOT_renormalized',greenlda,1,paw_dmft,pawprtvol=1,opt_wt=1) 166 167 !== Compare greenlda%occup and lda_occup: check that LDA green function is fine 168 !---------------------------------------------------------------------- 169 write(message,'(2a)') ch10,& 170 & ' == Check lda occ. mat. from green with respect to the direct calc ==' 171 call wrtout(std_out,message,'COLL') 172 call diff_oper("Occup from LDA green function",& 173 & "LDA occupations",greenlda%occup,lda_occup,1,paw_dmft%dmft_tolfreq) 174 ! write(message,'(2a)') ch10,& 175 !& ' ***** => Warning : diff_oper is suppressed for test' 176 ! call wrtout(std_out,message,'COLL') 177 write(message,'(2a)') ch10,& 178 & ' ***** => Calculation of Green function is thus correct without self ****' 179 call wrtout(std_out,message,'COLL') 180 call destroy_green(greenlda) 181 182 !== Orthonormalize psichi 183 !---------------------------------------------------------------------- 184 call timab(621,1,tsec) 185 opt_renorm=3 186 if(paw_dmft%nspinor==2.and.paw_dmft%dmft_solv==5) opt_renorm=2 ! necessary to use hybri_limit in qmc_prep_ctqmc 187 ! ought to be generalized in the future 188 if(paw_dmft%dmft_solv/=-1) then 189 call psichi_renormalization(cryst_struc,paw_dmft,pawang,opt=opt_renorm) 190 191 ! =========================================================================== 192 ! == re-construct LDA green function with new psichis 193 ! =========================================================================== 194 write(message,'(6a)')& 195 & ch10,' ==============================================================='& 196 & ,ch10,' ===== LDA Green Function Calculation with renormalized psichi'& 197 & ,ch10,' ==============================================================' 198 end if 199 call timab(621,2,tsec) 200 201 call wrtout(std_out,message,'COLL') 202 call icip_green("LDA_renormalized",cryst_struc,greenlda,& 203 & paw_dmft,pawang,pawprtvol,self) 204 !call print_green('LDA_renormalized',greenlda,1,paw_dmft,pawprtvol=1,opt_wt=1) 205 206 !Need to store idmftloop and set it to zero to avoid useless print_energy in ab_out 207 idmftloop=paw_dmft%idmftloop 208 paw_dmft%idmftloop=0 209 call compute_energy(cryst_struc,energies_dmft,greenlda,paw_dmft,pawprtvol,pawtab,self,occ_type=" lda",part='both') 210 paw_dmft%idmftloop=idmftloop 211 212 if(paw_dmft%dmft_prgn==1) then 213 if(paw_dmft%lpsichiortho==1) then 214 call local_ks_green(greenlda,paw_dmft,prtopt=1) 215 end if 216 end if 217 call destroy_self(self) 218 !call printocc_green(greenlda,9,paw_dmft,3,chtype="LDA GREEN PSICHI") 219 220 write(message,'(6a)')& 221 & ch10,' ==============================================================='& 222 & ,ch10,' ===== Define Interaction and self-energy' & 223 & ,ch10,' ==============================================================' 224 call wrtout(std_out,message,'COLL') 225 !== define Interaction from input upawu and jpawu 226 !---------------------------------------------------------------------- 227 ABI_DATATYPE_ALLOCATE(hu,(cryst_struc%ntypat)) 228 call init_hu(cryst_struc,pawtab,hu,paw_dmft%dmftqmc_t2g) 229 call initialize_self(self,paw_dmft) 230 231 ! Set Hu in density representation for calculation of entropy if needed... 232 do itypat=1,cryst_struc%ntypat 233 if ( hu(itypat)%lpawu /= -1 ) then 234 call data4entropyDMFT_setHu(paw_dmft%forentropyDMFT,itypat,hu(itypat)%udens(:,:)) 235 end if 236 end do 237 238 !== define self from scratch or file and double counting 239 !---------------------------------------------------------------------- 240 !- Self allocated 241 call dc_self(greenlda%charge_matlu,cryst_struc,hu,self,paw_dmft%dmft_dc,pawprtvol) 242 243 !- Read self or do self=hdc 244 if(paw_dmft%dmft_solv==4) then 245 ! write(std_out,*) "shift before rw_self",self%qmc_shift(1) 246 call make_qmcshift_self(cryst_struc,hu,self) 247 end if 248 call timab(627,1,tsec) 249 call rw_self(self,paw_dmft,prtopt=2,opt_rw=1,istep_iter=1000*istep) 250 call timab(627,2,tsec) 251 252 !== If QMC is used, and self energy is read for file, then 253 !== one does NOT shifts the self-energy because it was already shifted when writed, 254 !== and thus then weiss will be shifted 255 !---------------------------------------------------------------------- 256 !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf==1) & 257 !& call make_qmcshift_self(cryst_struc,hu,self) 258 !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf/=1) & 259 !& call make_qmcshift_self(cryst_struc,hu,self,apply=.true.) 260 261 call destroy_green(greenlda) ! destroy LDA green function 262 call print_self(self,"print_dc",paw_dmft,prtopt=2) 263 264 265 266 !=========================================================================== 267 !== construct green function with the self-energy. 268 !=========================================================================== 269 write(message,'(6a)') & 270 & ch10,' =================================================================', & 271 & ch10,' ===== Green Function Calculation with input self-energy ========', & 272 & ch10,' =================================================================' 273 call wrtout(std_out,message,'COLL') 274 call icip_green("Green_inputself",cryst_struc,green,& 275 & paw_dmft,pawang,pawprtvol,self,opt_self=1) 276 !call print_green('beforefermi_green',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 277 ! stop 278 ! call leave_new('COLL') 279 280 !== Find fermi level 281 !--------------------------------------------------------------------- 282 !write(message,'(2a,i3,13x,a)') ch10,' === Compute green function from self-energy' 283 call fermi_green(cryst_struc,green,paw_dmft,pawang,self) 284 !== define weiss field only for the local quantities (opt_oper=2) 285 !---------------------------------------------------------------------- 286 ! write(std_out,*) "nkpt befreo init_greenweiss",ifreq,paw_dmft%nkpt 287 call init_green(weiss,paw_dmft,opt_oper_ksloc=2) 288 ! do ifreq=1,weiss%nw 289 ! write(std_out,*) "nkpt from weiss1",ifreq,weiss%oper(ifreq)%nkpt 290 ! enddo 291 292 !== Check fourier transforms 293 !---------------------------------------------------------------------- 294 if(check==1) then 295 call check_fourier_green(cryst_struc,green,paw_dmft,pawang) 296 end if 297 298 !== If QMC is used, and self energy is not read for file, then 299 !== one shifts the self-energy, and thus then weiss will be shifted 300 !== after dyson, in a coherent way concerning qmc_shift and qmc_xmu. 301 !---------------------------------------------------------------------- 302 !if(paw_dmft%dmft_solv==4.and.paw_dmft%dmft_rslf/=1) & 303 !& call make_qmcshift_self(cryst_struc,hu,self,apply=.true.) 304 !if(paw_dmft%dmft_solv==4) write(std_out,*) "shift after make_qmcshift_self",self%qmc_shift(1) 305 306 write(message,'(6a)') & 307 & ch10,' ======================================================'& 308 & ,ch10,' ===== DMFT Loop starts here ========'& 309 & ,ch10,' ======================================================' 310 call wrtout(std_out,message,'COLL') 311 312 !======================================================================= 313 !=== dmft loop ======================================================= 314 do idmftloop=1, paw_dmft%dmft_iter 315 !paw_dmft%idmftloop=idmftloop 316 paw_dmft%idmftloop=paw_dmft%idmftloop+1 317 ! ======================================================================= 318 istep_iter=1000*istep+idmftloop 319 320 write(message,'(2a,i3,13x,a)') ch10,& 321 & ' ===== DMFT Loop : ITER number',paw_dmft%idmftloop,'========' 322 call wrtout(std_out,message,'COLL') 323 324 ! == Dyson Equation G,self -> weiss(w) 325 ! --------------------------------------------------------------------- 326 call dyson(green,paw_dmft,self,weiss,opt_weissself=1) 327 ! call print_green('afterDyson',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 328 ! stop 329 ! call leave_new('COLL') 330 331 ! == Printout local "occupations" from weiss field (useless) 332 if(abs(pawprtvol)>3) then 333 call integrate_green(cryst_struc,weiss,paw_dmft,& 334 & pawang,prtopt=2,opt_ksloc=2) 335 call printocc_green(weiss,5,paw_dmft,3,opt_weissgreen=1) 336 end if 337 338 ! === Prepare data, solve Impurity problem: weiss(w) -> G(w) 339 ! --------------------------------------------------------------------- 340 call initialize_self(self_new,paw_dmft) 341 342 call impurity_solve(cryst_struc,green,hu,& 343 & paw_dmft,pawang,pawtab,self,self_new,weiss,pawprtvol) ! weiss-> green, or self if dmft_solv=1 344 ! if(paw_dmft%dmft_solv==4) write(std_out,*) "shift after impurity",self%qmc_shift(1) 345 346 ! == Compute double counting from charge from green_solver 347 ! --------------------------------------------------------------------- 348 if (green%has_charge_matlu_solver/=2) green%charge_matlu_solver=green%charge_matlu 349 call dc_self(green%charge_matlu_solver,cryst_struc,hu,self_new,paw_dmft%dmft_dc,pawprtvol) 350 351 ! == Solve dyson equation. G_imp(w), weiss_imp(w) -> Self_imp(w) 352 ! --------------------------------------------------------------------- 353 ! if dmft_solv==1, self is computed previously 354 if(abs(paw_dmft%dmft_solv)/=1) then 355 call dyson(green,paw_dmft,self_new,weiss,opt_weissself=2) 356 end if 357 ! do ifreq=1,green%nw 358 ! call sym_matlu(cryst_struc,self%oper(ifreq)%matlu,pawang) 359 ! enddo 360 361 ! == Possibility if imposing self (opt_rw==3) 362 ! --------------------------------------------------------------------- 363 call timab(627,1,tsec) 364 call rw_self(self_new,paw_dmft,prtopt=2,opt_rw=3,istep_iter=istep_iter) 365 call timab(627,2,tsec) 366 367 ! print dc just computed before and self computed in before in dyson or 368 ! impurity_solve 369 if(abs(pawprtvol)>=3) then 370 write(message,'(2a)') ch10," == New self" 371 call wrtout(std_out,message,'COLL') 372 call print_self(self_new,"print_dc",paw_dmft,2) 373 write(message,'(2a)') ch10," == Old self" 374 call wrtout(std_out,message,'COLL') 375 call print_self(self,"print_dc",paw_dmft,2) 376 end if 377 378 ! if(paw_dmft%dmft_solv==4) write(std_out,*) "shift before computeenergy ",self%qmc_shift(1) 379 ! == Compute Energy with NEW self-energy and edc from green_solver, 380 ! new local green function and old occupations for eband 381 ! fermi level not optimized for this self_energy. 382 ! --------------------------------------------------------------------- 383 ! green= local green function and local charge comes directly from solver 384 ! green= ks green function and occupations comes from old_self 385 call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,& 386 & pawtab,self_new,occ_type="nlda",part=part2) 387 388 ! == Mix new and old self_energies and double countings 389 ! --------------------------------------------------------------------- 390 call new_self(self,self_new,paw_dmft,1) ! self,self_new => self 391 write(message,'(2a)') ch10," == After mixing," 392 !print *, " my_rank newself", my_rank,self%oper(1)%matlu(1)%mat(1,1,1,1,1) 393 call wrtout(std_out,message,'COLL') 394 call print_self(self,"print_dc",paw_dmft,2) ! print self and DC 395 call destroy_self(self_new) 396 397 ! == Compute green function self -> G(k) 398 ! --------------------------------------------------------------------- 399 call compute_green(cryst_struc,green,paw_dmft,pawang,1,self,opt_self=1,opt_nonxsum=1) 400 401 call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=3,opt_ksloc=3,opt_diff=1) !,opt_nonxsum=1) 402 403 call printocc_green(green,5,paw_dmft,3,chtype="DMFT") 404 ! call printocc_green(green,9,paw_dmft,3,chtype="DMFT FULL") 405 if(paw_dmft%lpsichiortho==1.and.paw_dmft%dmft_prgn==1) then 406 call local_ks_green(green,paw_dmft,prtopt=1) 407 end if 408 409 ! == Find fermi level 410 ! --------------------------------------------------------------------- 411 call fermi_green(cryst_struc,green,paw_dmft,pawang,self) 412 ! call leave_new('COLL') 413 414 ! == Compute Energy with Mixed self-energy and green function recomputed with new self 415 ! --------------------------------------------------------------------- 416 ! green= lattice green function computed from self for a given chemical potential mu (self_mixed,mu) 417 ! green= local green function is computed from lattice green function(self_mixed,mu) 418 ! green= occupations are computed with lattice green function(self_mixed,mu) 419 call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,occ_type="nlda",part=part3) 420 421 ! == Save self on disk 422 ! --------------------------------------------------------------------- 423 call timab(627,1,tsec) 424 call rw_self(self,paw_dmft,prtopt=2,opt_rw=2) 425 call timab(627,2,tsec) 426 427 ! == Test convergency 428 ! --------------------------------------------------------------------- 429 char_enddmft="DMFT (end of DMFT loop)" 430 if(green%ifermie_cv==1.and.self%iself_cv==1.and.green%ichargeloc_cv==1.and.paw_dmft%idmftloop>1) then 431 write(message,'(a,8x,a)') ch10,"DMFT Loop is converged !" 432 call wrtout(std_out,message,'COLL') 433 char_enddmft="converged DMFT" 434 exit 435 end if 436 ! ======================================================================= 437 ! === end dmft loop ==================================================== 438 end do 439 !========================================================================= 440 441 !== Save self on disk 442 !------------------------------------------------------------------------- 443 call timab(627,1,tsec) 444 call rw_self(self,paw_dmft,prtopt=2,opt_rw=2) 445 call timab(627,2,tsec) 446 447 !paw_dmft%idmftloop=0 448 449 write(message,'(2a,13x,a)') ch10,' ===== DMFT Loop : END ',& 450 & '========' 451 call wrtout(std_out,message,'COLL') 452 453 ! compute Edc for U=1 and J=U/J 454 call init_energy(cryst_struc,energies_tmp) 455 !call compute_ldau_energy(cryst_struc,energies_tmp,green,paw_dmft,pawtab) 456 call compute_ldau_energy(cryst_struc,energies_tmp,green,paw_dmft,pawtab,paw_dmft%forentropyDMFT%J_over_U) 457 call data4entropyDMFT_setDc(paw_dmft%forentropyDMFT,energies_tmp%e_dc(:)) 458 call destroy_energy(energies_tmp,paw_dmft) 459 460 !== Compute final values for green functions, occupations, and spectral function 461 !-------------------------------------------------------------------------------- 462 !Do not compute here, because, one want a energy computed after the 463 !solver (for Hubbard I and LDA+U). 464 call compute_green(cryst_struc,green,paw_dmft,pawang,1,self,opt_self=1,opt_nonxsum=1) 465 call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=2,opt_ksloc=3) !,opt_nonxsum=1) 466 !call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,opt=0) 467 idmftloop=paw_dmft%idmftloop 468 paw_dmft%idmftloop=0 469 call compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,occ_type="nlda",part="band") 470 paw_dmft%idmftloop=idmftloop 471 472 !write(message,'(2a,13x,a)') ch10,' ===== DMFT Loop is finished' 473 !call wrtout(ab_out,message,'COLL') 474 !write(std_out,*) "PRINTOCC INITIAL" 475 call printocc_green(green,9,paw_dmft,3,chtype=char_enddmft) 476 !write(std_out,*) "KS=czero" 477 !green%occup%ks=czero 478 !write(std_out,*) "PRINTOCC AFTER KS=0" 479 !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT") 480 !write(std_out,*) "UPFOLD_OPER" 481 !call upfold_oper(green%occup,paw_dmft,1) 482 !write(std_out,*) "PRINTOCC AFTER UPFOLD_OPER" 483 !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT") 484 !write(std_out,*) "MATLU=czero" 485 !green%occup%matlu(1)%mat=czero 486 !green%occup%ks(:,:,:,:)=cmplx(real(green%occup%ks(:,:,:,:))) 487 !write(std_out,*) "PRINTOCC AFTER MATLU=0 AND IMAG KS=0" 488 !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT") 489 !write(std_out,*) "LOC_OPER" 490 !call loc_oper(green%occup,paw_dmft,1) 491 !write(std_out,*) "PRINTOCC AFTER LOC_OPER" 492 !call printocc_green(green,9,paw_dmft,3,chtype="converged DMFT") 493 !call flush_unit(std_out) 494 !call leave_new('COLL') 495 if(paw_dmft%dmft_solv<=2.and.paw_dmft%prtdos>=1) then 496 call spectral_function(cryst_struc,green,hu,& 497 & paw_dmft,pawang,pawtab,self,pawprtvol) 498 end if 499 call destroy_green(weiss) 500 call destroy_green(green) 501 !todo_ab rotate back density matrix into unnormalized basis just for 502 !printout 503 call destroy_hu(hu,cryst_struc%ntypat,paw_dmft%dmftqmc_t2g) 504 call destroy_self(self) 505 call destroy_energy(energies_dmft,paw_dmft) 506 507 write(message,'(2a,13x,a)') ch10,' ===== DMFT : END ',& 508 & '========' 509 call wrtout(std_out,message,'COLL') 510 511 ABI_DATATYPE_DEALLOCATE(hu) 512 513 DBG_EXIT("COLL") 514 515 end subroutine dmft_solve