TABLE OF CONTENTS


ABINIT/dmft_solve [ Functions ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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