TABLE OF CONTENTS


ABINIT/dmft_solve [ Functions ]

[ Top ] [ 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
      m_spectral_function,timab,wrtout

SOURCE

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

ABINIT/dyson [ Functions ]

[ Top ] [ Functions ]

NAME

 dyson

FUNCTION

 Use the Dyson Equation to compute self-energy from green function

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

  dtset <type(dataset_type)>=all input variables for this dataset
  istep    =  step of iteration for LDA.
  lda_occup
  mpi_enreg=informations about MPI parallelization
  paw_dmft =  data for self-consistent LDA+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 LDA+DMFT calculations.

NOTES

PARENTS

      dmft_solve

CHILDREN

      add_matlu,copy_green,destroy_green,init_green,inverse_oper,timab,wrtout

SOURCE

891 subroutine dyson(green,paw_dmft,self,weiss,opt_weissself)
892 
893  use defs_basis
894  use m_abicore
895  use m_errors
896 
897  use m_time,         only : timab
898  use m_paw_dmft, only: paw_dmft_type
899  use m_crystal, only : crystal_t
900  use m_green, only : green_type, destroy_green,init_green,copy_green
901  use m_oper, only : oper_type,inverse_oper
902  use m_matlu, only : matlu_type,add_matlu,print_matlu
903  use m_self, only : self_type
904 
905 !This section has been created automatically by the script Abilint (TD).
906 !Do not modify the following lines by hand.
907 #undef ABI_FUNC
908 #define ABI_FUNC 'dyson'
909 !End of the abilint section
910 
911  implicit none
912 
913 !Arguments ------------------------------------
914 !scalars
915  type(green_type),intent(inout)  :: green
916  type(paw_dmft_type), intent(inout) :: paw_dmft
917  type(self_type),intent(inout)  :: self
918  type(green_type),intent(inout)  :: weiss
919  integer,intent(in) :: opt_weissself
920 ! type(paw_dmft_type), intent(inout)  :: paw_dmft
921 
922 !Local variables ------------------------------
923 !arrays
924  real(dp) :: tsec(2)
925 !scalars
926  integer :: ifreq,natom,nspinor,nsppol,weissinv
927  type(green_type)  :: greeninv
928  character(len=500) :: message
929 ! type
930 ! type(matlu_type), pointer :: matlutemp,matlu1,matlu2
931 !************************************************************************
932 
933  call timab(623,1,tsec)
934  DBG_ENTER("COLL")
935  natom=green%oper(1)%natom
936  nsppol=green%oper(1)%nsppol
937  nspinor=green%oper(1)%nspinor
938  weissinv=1
939  if(opt_weissself==1) then
940    write(message,'(2a,i3,13x,a)') ch10,'  ===  Use Dyson Equation => weiss '
941    call wrtout(std_out,message,'COLL')
942  else if(opt_weissself==2) then
943    write(message,'(2a,i3,13x,a)') ch10,'  ===  Use Dyson Equation => self'
944    call wrtout(std_out,message,'COLL')
945  end if
946 
947 !call xmpi_barrier(spaceComm)
948  if(paw_dmft%dmft_solv==2) weissinv=0
949  call init_green(greeninv,paw_dmft,opt_oper_ksloc=2,wtype=green%w_type)
950  call copy_green(green,greeninv,opt_tw=2)
951 
952  do ifreq=1,green%nw
953    if(opt_weissself==1) then
954      call inverse_oper(greeninv%oper(ifreq),option=1,prtopt=1)
955 !    warning green is now inversed
956      call add_matlu(greeninv%oper(ifreq)%matlu,self%oper(ifreq)%matlu,&
957 &     weiss%oper(ifreq)%matlu,natom,sign_matlu2=1)
958      call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1)
959    else if(opt_weissself==2) then
960 
961    ! 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))
962    ! 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))
963    ! 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))
964 !    call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1)
965 !    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))
966 !    call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1)
967 !    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))
968 
969 !    write(std_out,*) "-----------------------IFREQ",ifreq
970 !    call print_matlu(greeninv%oper(ifreq)%matlu,paw_dmft%natom,1,opt_diag=-1)
971      call inverse_oper(greeninv%oper(ifreq),option=1,prtopt=1)
972 !    call print_matlu(greeninv%oper(ifreq)%matlu,paw_dmft%natom,1,opt_diag=-1)
973 !    If paw_dmft%dmft_solv==2, then inverse of weiss function is
974 !    computed in m_hubbard_one.F90
975      if(weissinv/=0) then
976        call inverse_oper(weiss%oper(ifreq),option=1,prtopt=1)
977      end if
978 
979 !    write(std_out,*) weiss%oper(1)%matlu(ifreq)%mat(1,1,1,1,1),"-",greeninv%oper(ifreq)
980      call add_matlu(weiss%oper(ifreq)%matlu,greeninv%oper(ifreq)%matlu,&
981 &     self%oper(ifreq)%matlu,natom,sign_matlu2=-1)
982 
983    ! 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))
984    ! 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))
985    ! 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))
986    else
987      message = " BUG in dyson.F90"
988      MSG_BUG(message)
989    end if
990  end do
991 
992  call destroy_green(greeninv)
993 
994 
995  call timab(623,2,tsec)
996  DBG_EXIT("COLL")
997 
998 end subroutine dyson

ABINIT/impurity_solve [ Functions ]

[ Top ] [ Functions ]

NAME

 impurity_solve

FUNCTION

 Solve the Impurity problem

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
  lda_occup
  paw_dmft =  data for self-consistent LDA+DMFT calculations.
  pawang <type(pawang)>=paw angular mesh and related data
  pawtab <type(pawtab)>

OUTPUT

  paw_dmft =  data for self-consistent LDA+DMFT calculations.

NOTES

PARENTS

      m_dmft

CHILDREN

      copy_green,destroy_green_tau,flush_unit,fourier_green,m_hubbard_one
      init_green_tau,integrate_green,ldau_self,print_green,print_matlu
      printocc_green,qmc_prep_ctqmc,timab,trace_oper,wrtout

SOURCE

589 subroutine impurity_solve(cryst_struc,green,hu,paw_dmft,&
590 & pawang,pawtab,self_old,self_new,weiss,pawprtvol)
591 
592 
593  use defs_basis
594  use defs_abitypes
595  use m_errors
596  use m_abicore
597 
598  use m_time,    only : timab
599  use m_crystal, only : crystal_t
600  use m_green, only : green_type, fourier_green&
601 & ,init_green_tau,destroy_green_tau,print_green,printocc_green,integrate_green,copy_green
602  use m_paw_dmft, only : paw_dmft_type
603  use m_oper,     only : trace_oper
604  use m_hu, only : hu_type
605  use m_matlu, only : matlu_type,print_matlu,init_matlu,destroy_matlu
606  use m_self, only : self_type
607  use m_energy, only : energy_type
608  use m_pawang, only : pawang_type
609  use m_pawtab, only : pawtab_type
610  use m_io_tools, only : flush_unit
611  use m_hubbard_one, only : hubbard_one
612  use m_ldau_self, only : ldau_self
613  use m_forctqmc, only : qmc_prep_ctqmc
614 
615 !This section has been created automatically by the script Abilint (TD).
616 !Do not modify the following lines by hand.
617 #undef ABI_FUNC
618 #define ABI_FUNC 'impurity_solve'
619 !End of the abilint section
620 
621  implicit none
622 
623 !Arguments ------------------------------------
624 !scalars
625 ! type(pawang_type), intent(in) :: pawang
626  type(crystal_t),intent(in) :: cryst_struc
627  type(green_type), intent(inout) :: weiss
628  type(green_type), intent(inout) :: green !vz_i
629  type(hu_type),intent(inout) :: hu(cryst_struc%ntypat)
630  !type(MPI_type), intent(in) :: mpi_enreg
631  type(pawang_type), intent(in) :: pawang
632  type(pawtab_type),intent(in)  :: pawtab(cryst_struc%ntypat)
633  type(paw_dmft_type), intent(inout)  :: paw_dmft
634  type(self_type), intent(inout) :: self_new
635  type(self_type), intent(inout) :: self_old
636  integer, intent(in) :: pawprtvol
637 
638 !Local variables ------------------------------
639  real(dp) :: tsec(2)
640  character(len=500) :: message
641  complex(dpc) :: xx
642  integer :: ifreq
643 ! integer iatom,il,i_nd,isppol,lpawu,im,Nd,nrat,nsweeptot
644 ! real(dp) :: acc,kx
645 ! real(dp), allocatable :: correl(:,:),g0(:,:),gtmp(:,:)
646 !scalars
647 !************************************************************************
648 !character(len=500) :: message
649 
650  call timab(622,1,tsec)
651 !=======================================================================
652 !== Prepare data for Hirsch Fye QMC
653 !== NB: for CTQMC, Fourier Transformation are done inside the CTQMC code
654 !=======================================================================
655  if(abs(paw_dmft%dmft_solv)==4) then
656 !  == Initialize weiss and green functions for fourier transformation
657 !  -------------------------------------------------------------------
658    write(message,'(2a,i3,13x,a)') ch10,'   ===  Initialize Weiss field G_0(tau)'
659    call wrtout(std_out,message,'COLL')
660    call init_green_tau(weiss,paw_dmft)
661    call init_green_tau(green,paw_dmft)
662 !  in init_solver
663 
664 !  == Print weiss function G_0(tau=0-) before computation (really useless check)
665 !  ------------------------------------------------------------------------------
666    if(abs(pawprtvol)>3) then
667      write(message,'(2a,i3,13x,a)') ch10,'   ===  Check G_0(tau=0-) first'
668      call wrtout(std_out,message,'COLL')
669      call printocc_green(weiss,6,paw_dmft,3)
670    end if
671 
672 !  == Fourier transform of weiss Field
673 !  ------------------------------------
674 !  for fourier of KS green functions
675 !  call fourier_green(cryst_struc,weiss,mpi_enreg,paw_dmft,pawang,pawtab,1)
676    write(message,'(2a,i3,13x,a)') ch10,'   ===  Inverse Fourier Transform w->t of Weiss Field'
677    call wrtout(std_out,message,'COLL')
678    call fourier_green(cryst_struc,weiss,paw_dmft,pawang,opt_ksloc=2,opt_tw=-1)
679 
680 !  == Print weiss function G2_0(tau=0-)
681 !  --------------------------------------
682    call printocc_green(weiss,6,paw_dmft,3,opt_weissgreen=1)
683 
684 !  for fourier of KS green functions
685 !  call fourier_green(cryst_struc,weiss,mpi_enreg,paw_dmft,pawang,pawtab,1)
686 !  == Print G_0(tau) in files
687 !  ---------------------------
688    if(paw_dmft%dmft_prgn==1) then
689      call print_green('weiss',weiss,1,paw_dmft,pawprtvol=1,opt_wt=2)
690    end if
691 
692  else if(abs(paw_dmft%dmft_solv)>=5) then
693 !  == Initialize  green functions for imaginary times
694 !  -------------------------------------------------------------------
695    write(message,'(2a,i3,13x,a)') ch10,'   ===  Initialize Green function G(tau)'
696    call wrtout(std_out,message,'COLL')
697    call init_green_tau(green,paw_dmft)
698 
699  end if
700 !=======================================================================
701 !== End preparation of QMC
702 !=======================================================================
703 
704 !=======================================================================
705 !== Solve impurity model   =============================================
706 !=======================================================================
707  write(message,'(2a,i3,13x,a)') ch10,'  ===  Solve impurity model'
708  call wrtout(std_out,message,'COLL')
709  if(abs(paw_dmft%dmft_solv)==1) then
710 
711 !  == LDA+U for test -> self
712 !  -------------------
713    call ldau_self(cryst_struc,green,paw_dmft,&
714 &   pawtab,self_new,opt_ldau=1,prtopt=pawprtvol)
715  else if(abs(paw_dmft%dmft_solv)==2) then
716 
717 !  == Hubbard One -> green
718 !  -------------------
719    call hubbard_one(cryst_struc,green,hu,paw_dmft,&
720 &   pawang,pawprtvol,self_old%hdc,weiss)
721 
722  else if(abs(paw_dmft%dmft_solv)==4) then
723 
724 !  == QMC
725 !  -------------------
726    call copy_green(weiss,green,opt_tw=1)
727 !  call qmc_prep
728    message = '  ===  QMC not yet distributed '
729    MSG_ERROR(message)
730 !   call qmc_prep(cryst_struc,green,hu,mpi_enreg,paw_dmft,pawang&
731 !&   ,pawprtvol,self_old%qmc_xmu,weiss)
732 
733  else if(abs(paw_dmft%dmft_solv)>=5) then
734 
735 !  == Nothing
736 !  -------------------
737 !   call copy_green(weiss,green,opt_tw=1)
738 !   call copy_green(weiss,green,opt_tw=2)
739 
740    call qmc_prep_ctqmc(cryst_struc,green,self_old,hu,paw_dmft,pawang,pawprtvol,weiss)
741 
742 
743  else if(abs(paw_dmft%dmft_solv)==0) then
744 
745 !  == Nothing
746 !  -------------------
747 !  weiss%occup%has_operks=0 -> only local part is duplicated
748    call copy_green(weiss,green,opt_tw=2)
749  end if
750 !call print_green("invWeiss",cryst_struc,weiss,3,paw_dmft,pawtab,2)
751 
752 !=======================================================================
753 !== Treat data from HF QMC
754 !=======================================================================
755  if(abs(paw_dmft%dmft_solv)>=4) then
756 !  propagate qmc_shift (useful for compute_energy)
757    if(abs(paw_dmft%dmft_solv)==4) then
758      self_new%qmc_shift(:)=self_old%qmc_shift(:)
759      self_new%qmc_xmu(:)=self_old%qmc_xmu(:)
760    end if
761 
762 !  == Print local occupations from G(tau)
763 !  ---------------------------------------
764 
765 !  == Fourier back transform of green function G(tau)->G(iw_n) and
766 !  == compute occupations from g(tau)
767 !  -------------------------------------------------------------------
768    if(abs(paw_dmft%dmft_solv)==4) then
769      write(message,'(2a,i3,13x,a)') ch10,'   ===  Direct Fourier Transform t->w of Green Function'
770      call wrtout(std_out,message,'COLL')
771      call fourier_green(cryst_struc,green,paw_dmft,&
772 &     pawang,opt_ksloc=2,opt_tw=1)
773      do ifreq=1,green%nw
774        xx= green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
775        write(112,*) paw_dmft%omega_lo(ifreq),real(one/xx),aimag(one/xx)
776        write(113,*) paw_dmft%omega_lo(ifreq),real(xx),aimag(xx)
777      end do
778      call flush_unit(112)
779      call flush_unit(113)
780 !   if(paw_dmft%dmft_solv==5) stop
781      if(pawprtvol>=3) then
782        write(message,'(a,2x,a,f13.5)') ch10,&    ! debug
783 &      " == Print green function for small freq after fourier " ! debug
784        call wrtout(std_out,message,'COLL')    ! debug
785        call print_matlu(green%oper(1)%matlu,paw_dmft%natom,1)    ! debug
786      end if
787 
788      write(message,'(2a,i3,13x,a)') ch10,'   INVERSE FOURIER OF G0 SUPPRESSED'
789      call wrtout(std_out,message,'COLL')
790    end if
791    if(abs(paw_dmft%dmft_solv)==888) then
792 !  == Back fourier transform of G_0(tau) for compensation (try or comment or improve FT).
793 !  -------------------------------------------------------------------
794      write(message,'(2a,i3,13x,a)') ch10,'   ===  Direct Fourier transform t->w of Weiss'
795      call wrtout(std_out,message,'COLL')
796      call fourier_green(cryst_struc,weiss,paw_dmft,&
797 &     pawang,opt_ksloc=2,opt_tw=1)
798 
799      if(pawprtvol>=3) then
800        write(message,'(a,2x,a,f13.5)') ch10,&    ! debug
801 &      " == Print weiss function for small freq after fourier " ! debug
802        call wrtout(std_out,message,'COLL')    ! debug
803        call print_matlu(weiss%oper(1)%matlu,paw_dmft%natom,1)    ! debug
804      end if
805      call destroy_green_tau(weiss)
806    end if
807 
808 !  == Destroy tau part of green
809 !  -------------------------------------------------------------------
810    call trace_oper(green%occup_tau,green%charge_ks,green%charge_matlu_solver,2)
811    green%has_charge_matlu_solver=2
812    call destroy_green_tau(green)
813  end if
814 !=======================================================================
815 !== End Treat data for QMC
816 !=======================================================================
817 
818 !=======================================================================
819 !== Integrate green function and printout occupations
820 !=======================================================================
821 !For dmft_solv=-1,0,or 1 , the green function was not yet computed: it
822 !cannot be integrated
823 !=======================================================================
824  if(paw_dmft%dmft_solv>=2.and.green%w_type=="imag") then
825 !  ==  Integrate G(iw_n)
826 !  ---------------------
827    write(message,'(2a,i3,13x,a)') ch10,'   ===  Integrate local part of green function'
828    call wrtout(std_out,message,'COLL')
829    call integrate_green(cryst_struc,green,paw_dmft,&
830 &   pawang,prtopt=2,opt_ksloc=2,opt_after_solver=1)
831 
832 !  == Print local occupations from integration of G(iw_n)
833 !  --------------------------------------------------------
834    call printocc_green(green,5,paw_dmft,3)
835 
836 !  == Print G_loc(w)
837 !  --------------------------------------------------------
838    if(paw_dmft%dmft_prgn==1) then
839      call print_green('DMFT_IMPURITY',green,1,paw_dmft,pawprtvol=1,opt_wt=1)
840    end if
841  end if
842 !stop
843 
844  if(abs(pawprtvol)>0) then
845  end if
846 
847  call timab(622,2,tsec)
848 end subroutine impurity_solve

ABINIT/m_dmft [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dmft

FUNCTION

COPYRIGHT

 Copyright (C) 2006-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 .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 
29 #include "abi_common.h"
30 
31 MODULE m_dmft
32 
33  use defs_basis
34  implicit none
35 
36  private
37 
38  public :: dmft_solve
39  public :: impurity_solve
40  public :: dyson
41  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

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
  green  <type(green_type)>= green function data
  hu  <type(hu_type)>= datatype of type hu
  paw_dmft =  data for self-consistent LDA+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 LDA+DMFT calculations.

NOTES

PARENTS

      m_dmft

CHILDREN

      compute_green,copy_green,copy_matlu,dc_self,destroy_green,destroy_self
      dyson,hubbard_one,init_green,initialize_self,ldau_self,print_green
      rw_self,wrtout

SOURCE

1039 subroutine spectral_function(cryst_struc,green,hu,paw_dmft,&
1040 & pawang,pawtab,self_old,prtopt)
1041 
1042  use defs_basis
1043  use defs_abitypes
1044  use m_errors
1045  use m_abicore
1046 
1047  use m_crystal, only : crystal_t
1048  use m_green, only : init_green,green_type,print_green,copy_green,compute_green,destroy_green
1049  use m_matlu, only : copy_matlu
1050  use m_paw_dmft, only : paw_dmft_type
1051  use m_hu, only : hu_type
1052  use m_self, only : self_type,initialize_self,dc_self,destroy_self,rw_self
1053  use m_energy, only : energy_type
1054  use m_pawang, only : pawang_type
1055  use m_pawtab, only : pawtab_type
1056  use m_hubbard_one, only : hubbard_one
1057  use m_ldau_self, only : ldau_self
1058 
1059 !This section has been created automatically by the script Abilint (TD).
1060 !Do not modify the following lines by hand.
1061 #undef ABI_FUNC
1062 #define ABI_FUNC 'spectral_function'
1063 !End of the abilint section
1064 
1065  implicit none
1066 
1067 !Arguments ------------------------------------
1068 !scalars
1069 ! type(pawang_type), intent(in) :: pawang
1070  type(crystal_t),intent(in) :: cryst_struc
1071  type(green_type), intent(in) :: green
1072  type(hu_type),intent(inout) :: hu(cryst_struc%ntypat)
1073  !type(MPI_type), intent(inout) :: mpi_enreg
1074  type(pawang_type), intent(in) :: pawang
1075  type(pawtab_type),intent(in)  :: pawtab(cryst_struc%ntypat)
1076  type(self_type), intent(inout) :: self_old
1077  type(paw_dmft_type), intent(inout)  :: paw_dmft
1078  integer, intent(in) :: prtopt
1079 
1080 !Local variables ------------------------------
1081  character(len=500) :: message
1082  type(green_type) :: greenr
1083  type(green_type) :: weissr
1084  type(self_type) :: selfr
1085 !scalars
1086 !************************************************************************
1087 !character(len=500) :: message
1088 
1089 !   opt_oper_ksloc=3 to be able to compute spectral function.
1090  call init_green(greenr,paw_dmft,opt_oper_ksloc=3,wtype="real")
1091  call init_green(weissr,paw_dmft,wtype="real")
1092  call copy_matlu(green%occup%matlu,greenr%occup%matlu,paw_dmft%natom)
1093  call initialize_self(selfr,paw_dmft,wtype="real")
1094 !=======================================================================
1095 !== Solve impurity model with green function for real frequency
1096 !=======================================================================
1097  write(message,'(2a,i3,13x,a)') ch10,'  ===  Write Spectral function'
1098  call wrtout(std_out,message,'COLL')
1099  if(abs(paw_dmft%dmft_solv)==1) then
1100 
1101 !  == LDA+U for test
1102 !  -------------------
1103    call ldau_self(cryst_struc,greenr,paw_dmft,&
1104 &   pawtab,selfr,opt_ldau=1,prtopt=prtopt)
1105  else if(abs(paw_dmft%dmft_solv)==2) then
1106 
1107 !  == Hubbard One
1108 !  -------------------
1109    call hubbard_one(cryst_struc,greenr,hu,paw_dmft,&
1110 &   pawang,prtopt,self_old%hdc,weissr)
1111 
1112  else if(abs(paw_dmft%dmft_solv)==4) then
1113 
1114 !  == Nothing
1115 !  -------------------
1116    message = "spectral_function: This section of code is disabled!"
1117    MSG_ERROR(message)
1118    call copy_green(weissr,greenr,opt_tw=1)
1119 
1120  else if(abs(paw_dmft%dmft_solv)>=5) then
1121 
1122 !  == Nothing
1123 !  -------------------
1124    MSG_ERROR("Stopping before copy_green")
1125    call copy_green(weissr,greenr,opt_tw=1)
1126 
1127  else if(abs(paw_dmft%dmft_solv)==0) then
1128 
1129 !  == Nothing
1130 !  -------------------
1131 !  weiss%occup%has_operks=0 -> only local part is duplicated
1132    call copy_green(weissr,greenr,opt_tw=2)
1133  end if
1134 
1135 !=======================================================================
1136 !== Integrate green function and printout occupations
1137 !For dmft_solv=-1,0,or 1 , the green function was not computed: it
1138 !cannot be integrated
1139 !=======================================================================
1140  call dc_self(green%charge_matlu_solver,cryst_struc,hu,selfr,paw_dmft%dmft_dc,prtopt)
1141  if(abs(paw_dmft%dmft_solv)/=1.and.paw_dmft%dmft_solv/=0) then
1142    call dyson(greenr,paw_dmft,selfr,weissr,opt_weissself=2)
1143  end if
1144  call compute_green(cryst_struc,greenr,paw_dmft,pawang,1,selfr,opt_self=1)
1145  call print_green("realw",greenr,4,paw_dmft,pawprtvol=3)
1146  call rw_self(selfr,paw_dmft,prtopt=2,opt_rw=2)
1147 
1148  if(abs(prtopt)>0) then
1149  end if
1150  call destroy_self(selfr)
1151  call destroy_green(weissr)
1152  call destroy_green(greenr)
1153 
1154 end subroutine spectral_function