TABLE OF CONTENTS


ABINIT/m_energy [ Modules ]

[ Top ] [ Modules ]

NAME

  m_energy

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_energy
32 
33  use defs_basis
34  use m_errors
35  use m_abicore
36 
37  use m_pawtab, only : pawtab_type
38  use m_paw_correlations, only : pawuenergy
39 
40  use m_green, only : green_type,icip_green,destroy_green,compa_occup_ks
41  use m_self, only : self_type,initialize_self,destroy_self,print_self,new_self,make_qmcshift_self
42  use m_paw_dmft, only : paw_dmft_type
43 
44  use m_matlu, only : matlu_type,init_matlu,prod_matlu,diag_matlu,destroy_matlu,conjg_matlu,&
45 & ln_matlu,add_matlu,zero_matlu,shift_matlu,copy_matlu,trace_matlu,print_matlu
46  use m_oper, only : trace_oper
47  use m_crystal, only : crystal_t
48 
49  implicit none
50 
51  private
52 
53  public :: init_energy
54  public :: compute_energy
55  public :: compute_migdal_energy
56  public :: compute_ldau_energy
57  public :: destroy_energy
58  public :: print_energy
59  public :: compute_noninterentropy

m_energy/compute_band_energy [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 compute_band_energy

FUNCTION

INPUTS

  energies_dmft <type(energy_type)> = DMFT energy structure data
  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  occ_type=  character ("lda" or "nlda") for printing.
  fcalc_lda= if present, compute free energy instead of total energy.
  fcalc_lda= 2 do a free energy calculation with the dmft fermie level
           = 2/3 subtract fermie to eigenvalues and free energy
           =1/3 fermie_lda level and free energy
  ecalc_lda= 1 subtract fermie to eigenvalues

OUTPUT

SIDE EFFECTS

  energies_dmft <type(energy_type)> = DMFT energy structure data

PARENTS

      m_energy

CHILDREN

      wrtout

SOURCE

526 subroutine compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_lda,fcalc_lda,ecalc_dmft)
527 
528 
529 !This section has been created automatically by the script Abilint (TD).
530 !Do not modify the following lines by hand.
531 #undef ABI_FUNC
532 #define ABI_FUNC 'compute_band_energy'
533 !End of the abilint section
534 
535  implicit none
536 
537 !Arguments ------------------------------------
538 !type
539  type(energy_type),intent(inout) :: energies_dmft
540  type(green_type),intent(in) :: green
541  type(paw_dmft_type), intent(inout) :: paw_dmft
542  character(len=4), intent(in) :: occ_type
543  integer, intent(in), optional :: ecalc_lda
544  integer, intent(in), optional :: fcalc_lda
545  integer, intent(in), optional :: ecalc_dmft
546 ! integer :: prtopt
547 
548 !Local variables-------------------------------
549  integer :: ib,ikpt,isppol
550  real(dp) :: beta , fermie_used, totch,totch2,totch3
551  character(len=500) :: message
552 ! *********************************************************************
553  write(message,'(2a)') ch10,"  == Compute Band Energy terms for DMFT "
554  call wrtout(std_out,message,'COLL')
555  beta=one/paw_dmft%temp
556 
557 ! == Compute Band Energy
558 ! -----------------------------------------------------------------------
559  energies_dmft%eband_lda=zero
560  energies_dmft%eband_dmft=zero
561  totch=zero
562  totch2=zero
563  totch3=zero
564  do isppol=1,paw_dmft%nsppol
565    do ikpt=1,paw_dmft%nkpt
566      do ib=1,paw_dmft%mbandc
567        if(present(fcalc_lda)) then
568          if (fcalc_lda==1.or.fcalc_lda==3) fermie_used=paw_dmft%fermie_lda
569          if (fcalc_lda==2.or.fcalc_lda==4) fermie_used=paw_dmft%fermie ! only for B3 terms
570          if((paw_dmft%eigen_lda(isppol,ikpt,ib)-fermie_used).ge.zero) then
571            energies_dmft%eband_lda=energies_dmft%eband_lda - &
572 &             log ( one + exp ( - beta* (paw_dmft%eigen_lda(isppol,ikpt,ib)-fermie_used)))*paw_dmft%wtk(ikpt)
573          else
574            energies_dmft%eband_lda=energies_dmft%eband_lda  &
575 &             -(log ( one + exp ( beta* (paw_dmft%eigen_lda(isppol,ikpt,ib)-fermie_used))))*paw_dmft%wtk(ikpt) &
576 &              + beta*paw_dmft%eigen_lda(isppol,ikpt,ib)*paw_dmft%wtk(ikpt)
577            if(fcalc_lda==3.or.fcalc_lda==2) then
578                    ! subtract fermi level, (useful to directly count the number of electrons)
579              energies_dmft%eband_lda=energies_dmft%eband_lda  &
580 &                                 - beta*fermie_used*paw_dmft%wtk(ikpt)
581              totch=totch+paw_dmft%wtk(ikpt)
582            endif
583          endif
584        else ! usual calculation: total non interacting energy
585          fermie_used=paw_dmft%fermie_lda
586 !            write(std_out,*) "isppol,ikpt,ib",isppol,ikpt,ib
587 !            write(std_out,*) "paw_dmft%eigen_lda",paw_dmft%eigen_lda(isppol,ikpt,ib)
588 !            write(std_out,*) green%occup%ks(isppol,ikpt,ib,ib)
589 !            write(std_out,*) occup_fd(paw_dmft%eigen_lda(isppol,ikpt,ib),paw_dmft%fermie,paw_dmft%temp)
590          if(present(ecalc_lda)) then
591            if(ecalc_lda==1.or.ecalc_lda==3) fermie_used=paw_dmft%fermie_lda
592            if(ecalc_lda==2.or.ecalc_lda==4) fermie_used=paw_dmft%fermie ! only for B3 terms
593            if(ecalc_lda==3.or.ecalc_lda==2) then
594              energies_dmft%eband_lda=energies_dmft%eband_lda- &
595 &               occup_fd(paw_dmft%eigen_lda(isppol,ikpt,ib),fermie_used,paw_dmft%temp)*&
596 &               fermie_used*paw_dmft%wtk(ikpt)
597              totch2=totch2+paw_dmft%wtk(ikpt)*occup_fd(paw_dmft%eigen_lda(isppol,ikpt,ib),fermie_used,paw_dmft%temp)
598            endif
599          endif
600          energies_dmft%eband_lda=energies_dmft%eband_lda+ &
601 &           occup_fd(paw_dmft%eigen_lda(isppol,ikpt,ib),fermie_used,paw_dmft%temp)*&
602 &           paw_dmft%eigen_lda(isppol,ikpt,ib)*paw_dmft%wtk(ikpt)
603        endif
604        energies_dmft%eband_dmft=energies_dmft%eband_dmft+ &
605 &         green%occup%ks(isppol,ikpt,ib,ib)*&
606 &         paw_dmft%eigen_lda(isppol,ikpt,ib)*paw_dmft%wtk(ikpt)
607          totch3=totch3+paw_dmft%wtk(ikpt)*green%occup%ks(isppol,ikpt,ib,ib)
608          if(present(ecalc_dmft)) then
609            energies_dmft%eband_dmft=energies_dmft%eband_dmft- &
610 &              green%occup%ks(isppol,ikpt,ib,ib)*&
611 &              paw_dmft%fermie*paw_dmft%wtk(ikpt)
612          endif
613      enddo
614    enddo
615  enddo
616  if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1)  energies_dmft%eband_lda=two*energies_dmft%eband_lda
617  if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1)  energies_dmft%eband_dmft=two*energies_dmft%eband_dmft
618  if(present(fcalc_lda)) then
619    energies_dmft%eband_lda=energies_dmft%eband_lda/beta
620    if(fcalc_lda==3.or.fcalc_lda==2) write(std_out,*) "compute_band_energy totch",totch
621  endif
622  if(present(ecalc_lda)) then
623    if(ecalc_lda==3.or.ecalc_lda==2) write(std_out,*) "compute_band_energy totch2",totch2
624  endif
625 ! write(std_out,*) "compute_band_energy totch3",totch3
626 
627  if (occ_type==" lda") then
628    if(abs(energies_dmft%eband_lda-energies_dmft%eband_dmft)>tol5) then
629      write(message,'(5x,a,a,a,15x,a,f12.6,a,15x,a,5x,f12.5)')  "Warning !:"&
630 &     ,"Differences between band energy from LDA occupations",ch10&
631 &     ,"and LDA green function is:",energies_dmft%eband_lda-energies_dmft%eband_dmft,ch10&
632 &     ,"which is larger than",tol5
633      call wrtout(std_out,message,'COLL')
634      write(message,'(a)') &
635 &     "   Action: increase number of frequencies, or reduce the number of high energies_dmft bands"
636      call wrtout(std_out,message,'COLL')
637    else
638      write(message,'(a,a,a,10x,a,f12.6,a,10x,a,5x,f12.5)')  "          "&
639 &     ,"Differences between band energy from LDA occupations",ch10&
640 &     ,"and LDA green function is:",energies_dmft%eband_lda-energies_dmft%eband_dmft,ch10&
641 &     ,"which is smaller than",tol5
642      call wrtout(std_out,message,'COLL')
643    endif
644  endif
645 
646 
647 end subroutine compute_band_energy

m_energy/compute_energy [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 compute_energy

FUNCTION

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
  energies_dmft <type(energy_type)> = DMFT energy structure data
  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  pawprtvol
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  self  <type(self_type)>= self energy function data
  occ_type=  character ("lda" or "nlda") for printing.

OUTPUT

SIDE EFFECTS

 energies_dmft <type(energy_type)> = DMFT energy structure data

PARENTS

      m_dmft

CHILDREN

      wrtout

SOURCE

344 subroutine compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,occ_type,part)
345 
346 
347 !This section has been created automatically by the script Abilint (TD).
348 !Do not modify the following lines by hand.
349 #undef ABI_FUNC
350 #define ABI_FUNC 'compute_energy'
351 !End of the abilint section
352 
353  implicit none
354 
355 !Arguments ------------------------------------
356 !type
357  type(energy_type),intent(inout) :: energies_dmft
358  type(crystal_t),intent(in) :: cryst_struc
359  type(green_type),intent(in) :: green
360  type(paw_dmft_type), intent(inout) :: paw_dmft
361  type(pawtab_type),intent(in)  :: pawtab(cryst_struc%ntypat)
362  type(self_type), intent(inout) :: self
363  integer, intent(in) :: pawprtvol
364  character(len=4), intent(in) :: occ_type
365  character(len=4), intent(in) :: part
366 ! integer :: prtopt
367 
368 !Local variables-------------------------------
369  integer :: iatom,lpawu
370  integer :: natom,nspinor,nsppol
371  real(dp) :: beta
372  character(len=500) :: message
373  real(dp), allocatable :: e_hu_migdal(:)
374  real(dp) :: e_hu_migdal_tot
375 ! *********************************************************************
376  if(part=='both') then
377    write(message,'(2a)') ch10,"  == Compute LDA+DMFT energy terms "
378    call wrtout(std_out,message,'COLL')
379  else if(part=='band') then
380    write(message,'(2a)') ch10,"  == Compute LDA+DMFT energy terms : Band energy terms"
381    call wrtout(std_out,message,'COLL')
382  else if(part=='corr') then
383    write(message,'(2a)') ch10,"  == Compute LDA+DMFT energy terms : Correlation energy terms only"
384    call wrtout(std_out,message,'COLL')
385  else if(part=='none') then
386  endif
387 
388 ! Only imaginary frequencies here
389  if(green%w_type=="real".or.self%w_type=="real") then
390    message = 'compute_energy not implemented for real frequency'
391    MSG_BUG(message)
392  endif
393  natom=cryst_struc%natom
394  nsppol=paw_dmft%nsppol
395  nspinor=paw_dmft%nspinor
396  beta=one/paw_dmft%temp
397 
398  if(part=='band'.or.part=='both') then
399 ! == Compute Band Energy Alternative version: two steps
400 !                 == Compute Tr[ln G^{-1}] and -Tr[(Self-hdc)G_dmft]
401 ! -----------------------------------------------------------------------
402    if(part=='band') then ! ie if thdyn="fcalc" in m_dmft.F90
403 !     call compute_B3(cryst_struc,energies_dmft,eband2,green,mpi_enreg,paw_dmft,2,pawang,self,occ_type,0)
404 !     write(message,'(2a,f10.6)') ch10,"Compute Band energy test KS           ",eband2
405 !     call wrtout(std_out,message,'COLL')
406 !     call wrtout(ab_out,message,'COLL')
407 !
408 !     call compute_B3(cryst_struc,energies_dmft,eband2,green,mpi_enreg,paw_dmft,2,pawang,self,occ_type,1)
409 !     write(message,'(2a,f10.6)') ch10,"Compute Band energy test Self statique",eband2
410 !     call wrtout(std_out,message,'COLL')
411 !     call wrtout(ab_out,message,'COLL')
412 !
413 !! == Compute Band Energy (classical)
414 !! -----------------------------------------------------------------------
415 !     call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,fcalc_lda=3)
416 !     write(std_out,*) paw_dmft%fermie_lda,paw_dmft%fermie
417 !     write(message,'(2a,f10.6)') ch10,"Compute Band energy ref  free lda -ef ",energies_dmft%eband_lda
418 !     call wrtout(std_out,message,'COLL')
419      call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_lda=1)
420 !     write(message,'(2a,f10.6)') ch10,"Compute Band energy ref  -ef          ",energies_dmft%eband_dmft
421 !     call wrtout(std_out,message,'COLL')
422 !! call wrtout(ab_out,message,'COLL')
423 !     write(message,'(2a,f10.6)') ch10,"Compute Band energy ref   lda         ",energies_dmft%eband_lda
424 !     call wrtout(std_out,message,'COLL')
425 !! if(occ_type=="nlda") eband2=energies_dmft%eband_dmft
426    else
427      call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_lda=1)
428    endif
429 
430  endif
431  if(part=='corr'.or.part=='both') then
432 
433 ! == Compute Correlation energy from Migdal formula
434 ! -----------------------------------------------------------------------
435    ABI_ALLOCATE(e_hu_migdal,(cryst_struc%natom))
436    e_hu_migdal(:) = zero
437    call compute_migdal_energy(cryst_struc,e_hu_migdal,e_hu_migdal_tot,green,paw_dmft,pawprtvol,self)
438    energies_dmft%e_hu_mig(:)= e_hu_migdal(:)
439    energies_dmft%e_hu_mig_tot = e_hu_migdal_tot
440 ! write(std_out,*) "MIGDAL",e_hu_migdal_tot,e_hu_migdal
441    ABI_DEALLOCATE(e_hu_migdal)
442 
443 ! == Compute Correlation energy from QMC correlations.
444 ! -----------------------------------------------------------------------
445    energies_dmft%e_hu_qmc_tot = zero
446    do iatom=1,natom
447      lpawu=paw_dmft%lpawu(iatom)
448      if(lpawu/=-1) then
449        if(paw_dmft%dmft_solv==4.or.paw_dmft%dmft_solv==5.or.paw_dmft%dmft_solv==8) then
450          energies_dmft%e_hu_qmc(iatom) = green%ecorr_qmc(iatom)
451          energies_dmft%e_hu_qmc_tot = energies_dmft%e_hu_qmc_tot + green%ecorr_qmc(iatom)
452        endif
453      endif ! lpawu
454    enddo ! iatom
455 
456 ! == Compute LDA+U interaction energy
457 ! -----------------------------------------------------------------------
458    call compute_ldau_energy(cryst_struc,energies_dmft,green,paw_dmft,pawtab)
459    if(abs(paw_dmft%dmft_solv)<=1) then
460      energies_dmft%e_hu= energies_dmft%e_hu_ldau
461      energies_dmft%e_hu_tot= energies_dmft%e_hu_ldau_tot
462      if((abs(energies_dmft%e_hu_tot-energies_dmft%e_hu_mig_tot).ge.0.000001).and.(occ_type/=" lda")) then
463        write(message,'(2a,2e18.8,a)') ch10,'   BUG: Migdal energy and LDA+U energy do not coincide',&
464 &       energies_dmft%e_hu_tot,energies_dmft%e_hu_mig_tot,occ_type
465        MSG_ERROR(message)
466      endif
467    else if(paw_dmft%dmft_solv==2.or.paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7) then
468      energies_dmft%e_hu= energies_dmft%e_hu_mig
469      energies_dmft%e_hu_tot= energies_dmft%e_hu_mig_tot
470      energies_dmft%e_hu_qmc_tot = energies_dmft%e_hu_tot
471    else if(paw_dmft%dmft_solv==4.or.paw_dmft%dmft_solv==5.or.paw_dmft%dmft_solv==8) then
472      if(paw_dmft%dmft_solv==8) then
473        write(message,'(2a)') ch10,"Warning, energy is recently computed, not checked"
474        call wrtout(std_out,message,'COLL')
475      endif
476      energies_dmft%e_hu= energies_dmft%e_hu_qmc
477      energies_dmft%e_hu_tot= energies_dmft%e_hu_qmc_tot
478    endif
479 !   energies_dmft%edmft=energies_dmft%e_hu_mig_tot-energies_dmft%e_dc_tot
480    energies_dmft%edmft=energies_dmft%e_hu_tot-energies_dmft%e_dc_tot
481 
482 
483  endif ! part
484 
485 ! if(part='corr'.or.part='both') then
486  if(part/='none') then
487    call print_energy(cryst_struc,energies_dmft,pawprtvol,pawtab,paw_dmft%idmftloop)
488  endif
489 ! write(message,'(2a)') ch10," == The LDA+U self-energy is == "
490 ! call wrtout(std_out,message,'COLL')
491 ! call print_oper(self%oper(1),5,paw_dmft,2)
492 ! a voir: energies_dmft%e_hu_tot = energies_dmft%e_hu_ldau_tot
493 
494 end subroutine compute_energy

m_energy/compute_ldau_energy [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 compute_ldau_energy

FUNCTION

  Initialize noccmmp from green%occup and compute LDA+U energy with it

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  renorm = if present change U->1 and J-> renorm just for pawuenergy
           renorm = J/U for the real values (does not depend on "lambda" entropy)

OUTPUT

PARENTS

      m_dmft,m_energy

CHILDREN

      pawpuenergy,wrtout

SOURCE

831 subroutine compute_ldau_energy(cryst_struc,energies_dmft,green,paw_dmft,pawtab,renorm)
832 
833 
834 !This section has been created automatically by the script Abilint (TD).
835 !Do not modify the following lines by hand.
836 #undef ABI_FUNC
837 #define ABI_FUNC 'compute_ldau_energy'
838 !End of the abilint section
839 
840  implicit none
841 
842 !Arguments ------------------------------------
843 !type
844  type(energy_type),intent(inout) :: energies_dmft
845  type(crystal_t),intent(in) :: cryst_struc
846  type(green_type),intent(in) :: green
847  type(paw_dmft_type), intent(in) :: paw_dmft
848  type(pawtab_type),target,intent(in)  :: pawtab(cryst_struc%ntypat)
849  real(dp), optional, intent(in) :: renorm(:)
850 ! integer :: prtopt
851 
852 !Local variables-------------------------------
853  integer :: iatom,idijeff,im,im1,ispinor,ispinor1,isppol,ldim,lpawu
854  integer :: nocc,nsploop,prt_pawuenergy
855  real(dp) :: upawu,jpawu
856  real(dp) :: eldaumdcdc,eldaumdc,e_ee,e_dc,e_dcdc,xe1,xe2
857  character(len=500) :: message
858 ! arrays
859  integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
860  real(dp),allocatable :: noccmmp(:,:,:,:),nocctot(:)
861  type(pawtab_type),pointer :: pawtab_
862 
863 ! *********************************************************************
864 
865 ! - allocations
866 ! -----------------------------------------------------------------------
867 
868  nsploop=max(paw_dmft%nsppol,paw_dmft%nspinor**2)
869  nocc=nsploop
870  e_ee=zero
871  e_dc=zero
872  e_dcdc=zero
873  eldaumdc=zero
874  eldaumdcdc=zero
875  isppol=0
876  ispinor=0
877  ispinor1=0
878 
879 ! - Loop and call to pawuenergy
880 ! -----------------------------------------------------------------------
881  do iatom=1,cryst_struc%natom
882    pawtab_ => pawtab(cryst_struc%typat(iatom))
883    lpawu=paw_dmft%lpawu(iatom)
884    if(lpawu.ne.-1) then
885      ldim=2*lpawu+1
886      
887      ABI_ALLOCATE(noccmmp,(2,2*pawtab_%lpawu+1,2*pawtab_%lpawu+1,nocc))
888      ABI_ALLOCATE(nocctot,(nocc))
889      noccmmp(:,:,:,:)=zero ; nocctot(:)=zero
890 
891 ! - Setup nocctot and noccmmp
892 ! -----------------------------------------------------------------------
893      nocctot(:)=zero ! contains nmmp in the n m representation
894 ! Begin loop over spin/spinors to initialize noccmmp
895      do idijeff=1,nsploop
896        if(nsploop==2) then
897          isppol=spinor_idxs(1,idijeff)
898          ispinor=1
899          ispinor1=1
900        else if(nsploop==4) then
901          isppol=1
902          ispinor=spinor_idxs(1,idijeff)
903          ispinor1=spinor_idxs(2,idijeff)
904        else if(nsploop==1) then
905          isppol=1
906          ispinor=1
907          ispinor1=1
908        else
909          write(message,'(2a)') " BUG in m_energy: nsploop should be equal to 1, 2 or 4"
910          call wrtout(std_out,message,'COLL')
911        endif
912 ! Initialize noccmmp
913        do im1 = 1 , ldim
914          do im = 1 ,  ldim
915             noccmmp(1,im,im1,idijeff)=real(green%occup%matlu(iatom)%mat(im1,im,isppol,ispinor,ispinor1))
916             noccmmp(2,im,im1,idijeff)=aimag(green%occup%matlu(iatom)%mat(im1,im,isppol,ispinor,ispinor1))
917 !            noccmmp(1,im,im1,idijeff)=real(green%occup%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
918 !            noccmmp(2,im,im1,idijeff)=imag(green%occup%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
919          enddo
920        enddo
921 ! Compute nocctot
922        if(green%has_charge_matlu_solver/=2) then
923          do im1=1,ldim
924            if(nsploop==4) then
925              if(idijeff<=2) then
926                nocctot(1)=nocctot(1)+&
927 &                noccmmp(1,im1,im1,idijeff)
928              endif
929            else if (nsploop==2.or.nsploop==1) then
930              nocctot(idijeff)=nocctot(idijeff)+&
931 &              noccmmp(1,im1,im1,idijeff)
932            end if
933          enddo
934        else
935          if(nsploop==4) then
936            nocctot(1)=green%charge_matlu_solver(iatom,2) !  total nb of elec for nspinor=2 is (iatom,2) !!
937            nocctot(2)=zero
938            nocctot(3)=zero
939            nocctot(4)=zero
940           else if (nsploop==2) then
941            nocctot(1)=green%charge_matlu_solver(iatom,1) !  first spin
942            nocctot(2)=green%charge_matlu_solver(iatom,2) !  second one
943           else if (nsploop==1) then
944            nocctot(1)=green%charge_matlu_solver(iatom,1)
945          end if
946        endif
947      enddo
948 
949      xe1=e_dc
950      xe2=e_ee
951     ! write(std_out,*)" nocctot(1)",nocctot(1),green%charge_matlu_solver(iatom,1)
952      eldaumdc = zero
953      eldaumdcdc = zero
954      if ( present(renorm) ) then
955        upawu = one
956        jpawu = renorm(iatom)
957        prt_pawuenergy=0
958      else
959        upawu = pawtab_%upawu
960        jpawu = pawtab_%jpawu
961        prt_pawuenergy=3
962      end if
963 
964      call pawuenergy(iatom,eldaumdc,eldaumdcdc,noccmmp,nocctot,prt_pawuenergy,pawtab_,&
965 &                    dmft_dc=paw_dmft%dmft_dc,e_ee=e_ee,e_dc=e_dc,e_dcdc=e_dcdc,&
966 &                    u_dmft=upawu,j_dmft=jpawu)
967 
968      energies_dmft%e_dc(iatom)=e_dc-xe1
969      energies_dmft%e_hu_ldau(iatom)=e_ee-xe2
970 
971      ABI_DEALLOCATE(noccmmp)
972      ABI_DEALLOCATE(nocctot)
973    endif ! lpawu/=-1
974  enddo
975 
976 ! - gather results
977 ! -----------------------------------------------------------------------
978  energies_dmft%e_dc_tot=e_dc ! todo_ab: here or not ?
979  energies_dmft%e_hu_ldau_tot=e_ee
980 
981 end subroutine compute_ldau_energy

m_energy/compute_migdal_energy [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 compute_migdal_energy

FUNCTION

INPUTS

  cryst_struc <type(crystal_t)> = crystal structure data
  energies_dmft <type(energy_type)> = DMFT energy structure data
  green  <type(green_type)>= green function data
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  pawprtvol
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  self  <type(self_type)>= self energy function data

OUTPUT

  e_hu_mig(natom)= Migdal energy for each atom.
  e_hu_mig_tot= Total Migdal energy.

PARENTS

      m_energy

CHILDREN

      wrtout

SOURCE

677 subroutine compute_migdal_energy(cryst_struc,e_hu_migdal,e_hu_migdal_tot,green,paw_dmft,pawprtvol,self)
678 
679 #ifdef FC_INTEL
680 !DEC$ NOOPTIMIZE
681 #endif
682 
683 !This section has been created automatically by the script Abilint (TD).
684 !Do not modify the following lines by hand.
685 #undef ABI_FUNC
686 #define ABI_FUNC 'compute_migdal_energy'
687 !End of the abilint section
688 
689  implicit none
690 
691 !Arguments ------------------------------------
692 !type
693  type(crystal_t),intent(in) :: cryst_struc
694  type(green_type),intent(in) :: green
695  type(paw_dmft_type), intent(in) :: paw_dmft
696  real(dp),intent(out) :: e_hu_migdal_tot
697  real(dp),intent(out) :: e_hu_migdal(paw_dmft%natom)
698  type(self_type), intent(in) :: self
699  integer, intent(in) :: pawprtvol
700 ! integer :: prtopt
701 
702 !Local variables-------------------------------
703  integer :: iatom,ifreq,im,im1,ispinor,ispinor1,isppol,lpawu
704  integer :: natom,ndim,nspinor,nsppol,nwlo
705  real(dp) :: beta
706  complex(dpc) :: xmig_1,xmig_2,xmig_3,se,shift
707  character(len=500) :: message
708 ! *********************************************************************
709 
710 ! Only imaginary frequencies here
711  if(green%w_type=="real".or.self%w_type=="real") then
712    message = 'compute_migdal_energy not implemented for real frequency'
713    MSG_BUG(message)
714  endif
715 
716 ! == Compute Correlation energy from Migdal formula
717 ! -----------------------------------------------------------------------
718  natom=cryst_struc%natom
719  nsppol=paw_dmft%nsppol
720  nspinor=paw_dmft%nspinor
721  beta=one/paw_dmft%temp
722  nwlo=green%nw
723  if (green%nw/=self%nw) then
724    message = 'self and green do not contains the same number of frequencies'
725    MSG_BUG(message)
726  endif
727 ! write(std_out,*) "beta",beta
728 
729  xmig_1=zero
730  xmig_2=zero
731  xmig_3=zero
732 
733  e_hu_migdal_tot = zero
734  do iatom=1,natom
735    shift=czero
736    if(paw_dmft%dmft_solv==4) shift=self%qmc_shift(iatom)+self%qmc_xmu(iatom)
737 !   write(std_out,*) "shiftttt",shift
738    lpawu=paw_dmft%lpawu(iatom)
739    if(lpawu/=-1) then
740      xmig_1=czero
741      xmig_2=czero
742      xmig_3=czero
743      ndim=2*lpawu+1
744      do isppol=1,nsppol
745        do ispinor = 1 , nspinor
746          do ispinor1 = 1, nspinor
747            do im=1,ndim
748              do im1=1,ndim
749                do ifreq=1,nwlo
750 !                write(std_out,*) ifreq,xmig_1,imag(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)),&
751 !&                  green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor )
752                  xmig_1=xmig_1 + j_dpc/beta*       &
753 &                aimag(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1))* &
754 &                      green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor )* &
755 &                      paw_dmft%wgt_wlo(ifreq)
756 !                 if(ispinor==ispinor1.and.im==im1) then
757                    se=(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)-  &
758 &                      self%oper (nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1))
759 !                 else
760 !                   se=self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)
761 !                 endif
762                  xmig_2=xmig_2 + one/beta*real(se)* &
763 &                      green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor )* &
764 &                      paw_dmft%wgt_wlo(ifreq)
765 !                 if(ispinor==ispinor1.nd.im==im1.and.ifreq==1) then
766                  if(ifreq==1) then
767                    xmig_3=xmig_3 + &
768 &                   real(self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)+shift)* &
769 &                         green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two
770 !                   write(std_out,*) "xmig_3",xmig_3
771 !                   write(std_out,*) "self",self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)
772 !                   write(std_out,*) "shift",shift
773 !                   write(std_out,*) "occup", green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two
774                  endif
775                enddo
776 !               if(ispinor==ispinor1.and.im==im1) then
777 !                 xmig_3=xmig_3 + &
778 !&                 real(self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1))* &
779 !!&                         green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two
780 !               endif
781              enddo
782            enddo
783          enddo
784        enddo
785      enddo
786      if(nsppol==1.and.nspinor==1) then
787        e_hu_migdal(iatom)=two*real(xmig_1+xmig_2+xmig_3)
788      else
789        e_hu_migdal(iatom)=real(xmig_1+xmig_2+xmig_3)
790      endif
791      e_hu_migdal_tot = e_hu_migdal_tot + e_hu_migdal(iatom)
792      if(abs(pawprtvol)>=3) then
793        write(message,'(2a,3(a,5x,a,2f12.6))')ch10,&
794 &         "  Interaction energy: Decomposition of Migdal energy",ch10,&
795 &         "xmig_1=",xmig_1,ch10,&
796 &         "xmig_2=",xmig_2,ch10,&
797 &         "xmig_3=",xmig_3
798        call wrtout(std_out,message,'COLL')
799      endif
800    endif ! lpawu
801  enddo
802 
803 end subroutine compute_migdal_energy

m_energy/compute_noninterentropy [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 compute_noninterentropy

FUNCTION

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
  green  <type(green_type)>= green function data  only for Tr(G(self-hdc))
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data

OUTPUT

SIDE EFFECTS

PARENTS

CHILDREN

      wrtout

SOURCE

1006 subroutine compute_noninterentropy(cryst_struc,green,paw_dmft)
1007 
1008 
1009 !This section has been created automatically by the script Abilint (TD).
1010 !Do not modify the following lines by hand.
1011 #undef ABI_FUNC
1012 #define ABI_FUNC 'compute_noninterentropy'
1013 !End of the abilint section
1014 
1015  implicit none
1016 
1017 !Arguments ------------------------------------
1018 !type
1019  type(crystal_t),intent(in) :: cryst_struc
1020  type(green_type),intent(in) :: green
1021  type(paw_dmft_type), intent(inout) :: paw_dmft
1022 
1023 !Local variables-------------------------------
1024  integer :: ib,ikpt,isppol,natom,nspinor,nsppol
1025  real(dp) :: beta,eig,fermi,s_1,s_2,occ1,occ2,f_1,e_1,f_1a,s_1a,e_2
1026  character(len=800) :: message
1027 ! *********************************************************************
1028  write(message,'(2a,i6)') ch10,"  == Compute T*Entropy for fermi level and DFT-KS-LDA eigenvalues "
1029  call wrtout(std_out,message,'COLL')
1030 
1031  natom=cryst_struc%natom
1032  nsppol=paw_dmft%nsppol
1033  nspinor=paw_dmft%nspinor
1034  beta=one/paw_dmft%temp
1035  s_1=zero
1036  s_1a=zero
1037  f_1=zero
1038  f_1a=zero
1039  e_1=zero
1040  e_2=zero
1041  s_2=zero
1042  do isppol=1,paw_dmft%nsppol
1043    do ikpt=1,paw_dmft%nkpt
1044      do ib=1,paw_dmft%mbandc
1045        eig=paw_dmft%eigen_lda(isppol,ikpt,ib)
1046        fermi=paw_dmft%fermie_lda
1047        fermi=paw_dmft%fermie
1048        occ1=occup_fd(eig,fermi,paw_dmft%temp)
1049        occ2=green%occup%ks(isppol,ikpt,ib,ib)
1050 !       write(std_out,*) occ1,occ2
1051 
1052 !        entropy from Fermi Dirac
1053        if((occ1.ge.tol9).and.((one-occ1).ge.tol9)) then
1054          s_1=s_1+(occ1*log(occ1)+(one-occ1)*log(one-occ1))*paw_dmft%wtk(ikpt)
1055 !       write(std_out,*) occ1,one-occ1,"p1",(occ1*log(occ1)+(one-occ1)*log(one-occ1))*paw_dmft%wtk(ikpt)
1056        endif
1057 
1058 !        Free energy from Fermi Dirac
1059        if((eig-fermi).ge.zero) then ! occ1 -> 0 ; 1-occ1 -> 1
1060          f_1=f_1-paw_dmft%wtk(ikpt)/beta*log(one+exp(-beta*(eig-fermi)))
1061          f_1a=f_1a+paw_dmft%wtk(ikpt)/beta*log(one-occ1)
1062          s_1a=s_1a+((one-occ1)*log(one-occ1)+occ1*(-beta*(eig-fermi)+log(one-occ1)))*paw_dmft%wtk(ikpt)
1063        else ! occ1  -> 1 , 1-occ1 -> 0
1064          f_1=f_1-paw_dmft%wtk(ikpt)/beta*(log(one+exp(beta*(eig-fermi)))-beta*(eig-fermi))
1065          f_1a=f_1a-paw_dmft%wtk(ikpt)/beta*(-log(occ1)-beta*(eig-fermi))
1066          s_1a=s_1a+(occ1*log(occ1)+(one-occ1)*(beta*(eig-fermi)+log(occ1)))*paw_dmft%wtk(ikpt)
1067        endif
1068 
1069 !        Internal energy from Fermi Dirac
1070        e_1=e_1+(eig-fermi)*paw_dmft%wtk(ikpt)*occ1
1071        e_2=e_2+(eig-fermi)*paw_dmft%wtk(ikpt)*occ2
1072 
1073 !        entropy from green function occupations.
1074        if((occ2.ge.tol9).and.((one-occ2).ge.tol9)) then
1075 !       write(std_out,*) occ2,one-occ2,"p2",(occ2*log(occ2)+(one-occ2)*log(one-occ2))*paw_dmft%wtk(ikpt)
1076          s_2=s_2+(occ2*log(occ2)+(one-occ2)*log(one-occ2))*paw_dmft%wtk(ikpt)
1077        endif
1078      enddo
1079    enddo
1080  enddo
1081  s_1=-s_1*paw_dmft%temp
1082  s_1a=-s_1a*paw_dmft%temp
1083  s_2=-s_2*paw_dmft%temp
1084 
1085  write(message,'(8(2a,e20.9))') &
1086 & ch10," T*Entropy from Fermi Dirac occupations    ", s_1,&
1087 & ch10," T*Entropy from Fermi Dirac occupations  2 ", s_1a,&
1088 & ch10," T*Entropy from Green function occupations ", s_2,&
1089 & ch10," Free energy      F                        ", f_1,&
1090 & ch10," Free energy      Fa                       ", f_1a,&
1091 & ch10," internal energy  U                        ", e_1,&
1092 & ch10," internal energy  U  from Gr Func Occ      ", e_2,&
1093 & ch10," U-F                                       ", e_1-f_1
1094  call wrtout(std_out,message,'COLL')
1095 
1096 
1097 end subroutine compute_noninterentropy

m_energy/destroy_energy [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 destroy_energy

FUNCTION

  deallocate energies_dmft

INPUTS

  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data

OUTPUT

PARENTS

      m_dmft

CHILDREN

      wrtout

SOURCE

193 subroutine destroy_energy(energies_dmft,paw_dmft)
194 
195 
196 !This section has been created automatically by the script Abilint (TD).
197 !Do not modify the following lines by hand.
198 #undef ABI_FUNC
199 #define ABI_FUNC 'destroy_energy'
200 !End of the abilint section
201 
202  implicit none
203 
204 !Arguments ------------------------------------
205 !scalars
206  type(energy_type),intent(inout) :: energies_dmft
207  type(paw_dmft_type), intent(inout) :: paw_dmft
208 !Local variables-------------------------------
209 ! *********************************************************************
210   paw_dmft%edmft=energies_dmft%edmft
211  if ( allocated(energies_dmft%e_dc) )   then
212    ABI_DEALLOCATE(energies_dmft%e_dc)
213  end if
214  if ( allocated(energies_dmft%e_hu) )   then
215    ABI_DEALLOCATE(energies_dmft%e_hu)
216  end if
217  if ( allocated(energies_dmft%e_hu_ldau) )  then
218    ABI_DEALLOCATE(energies_dmft%e_hu_ldau)
219  end if
220  if ( allocated(energies_dmft%e_hu_mig) )  then
221    ABI_DEALLOCATE(energies_dmft%e_hu_mig)
222  end if
223   if ( allocated(energies_dmft%e_hu_qmc) )  then
224    ABI_DEALLOCATE(energies_dmft%e_hu_qmc)
225  end if
226 
227 
228 end subroutine destroy_energy

m_energy/energy_type [ Types ]

[ Top ] [ m_energy ] [ Types ]

NAME

  energy_type

FUNCTION

  This structured datatype contains interaction matrices for the correlated subspace

SOURCE

 72  type, public :: energy_type ! for each typat
 73 
 74   real(dp) :: eband_lda
 75 
 76   real(dp) :: eband_dmft
 77 
 78   real(dp) :: e_dc_tot
 79 
 80   real(dp) :: e_hu_tot
 81 
 82   real(dp) :: e_hu_ldau_tot
 83 
 84   real(dp) :: e_hu_mig_tot
 85 
 86   real(dp) :: e_hu_qmc_tot
 87 
 88   real(dp) :: edmft
 89 
 90   real(dp) :: natom
 91 
 92   real(dp), allocatable :: e_dc(:)
 93 
 94   real(dp), allocatable :: e_hu(:)
 95 
 96   real(dp), allocatable :: e_hu_ldau(:)
 97 
 98   real(dp), allocatable :: e_hu_mig(:)
 99 
100   real(dp), allocatable :: e_hu_qmc(:)
101 
102  end type energy_type

m_energy/init_energy [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 init_energy

FUNCTION

  Allocate variables used in type energy_type.

INPUTS

 OUTPUTS
 energies_dmft  = structure of data for dmft of type energy_type

PARENTS

      m_dmft

CHILDREN

      wrtout

SOURCE

132 subroutine init_energy(cryst_struc,energies_dmft)
133 
134 
135 !This section has been created automatically by the script Abilint (TD).
136 !Do not modify the following lines by hand.
137 #undef ABI_FUNC
138 #define ABI_FUNC 'init_energy'
139 !End of the abilint section
140 
141  implicit none
142 
143 !Arguments ------------------------------------
144 !type
145  type(crystal_t),intent(in) :: cryst_struc
146  type(energy_type), intent(inout) :: energies_dmft
147 !Local variables ------------------------------------
148 !************************************************************************
149 
150  ABI_ALLOCATE(energies_dmft%e_dc,(cryst_struc%natom))
151  ABI_ALLOCATE(energies_dmft%e_hu,(cryst_struc%natom))
152  ABI_ALLOCATE(energies_dmft%e_hu_ldau,(cryst_struc%natom))
153  ABI_ALLOCATE(energies_dmft%e_hu_mig,(cryst_struc%natom))
154  ABI_ALLOCATE(energies_dmft%e_hu_qmc,(cryst_struc%natom))
155  energies_dmft%e_dc=zero
156  energies_dmft%e_hu=zero
157  energies_dmft%e_hu_ldau=zero
158  energies_dmft%e_hu_mig=zero
159  energies_dmft%e_hu_qmc=zero
160  energies_dmft%eband_lda=zero
161  energies_dmft%eband_dmft=zero
162  energies_dmft%e_dc_tot=zero
163  energies_dmft%e_hu_tot=zero
164  energies_dmft%e_hu_ldau_tot=zero
165  energies_dmft%e_hu_mig_tot=zero
166  energies_dmft%e_hu_qmc_tot=zero
167  energies_dmft%edmft=zero
168  energies_dmft%natom=cryst_struc%natom
169 
170 end subroutine init_energy

m_energy/occup_fd [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 occup_fd

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

      wrtout

SOURCE

1118  function occup_fd(eig,fermie,temp)
1119 
1120 
1121 !This section has been created automatically by the script Abilint (TD).
1122 !Do not modify the following lines by hand.
1123 #undef ABI_FUNC
1124 #define ABI_FUNC 'occup_fd'
1125 !End of the abilint section
1126 
1127  implicit none
1128 
1129 !Arguments ------------------------------------
1130 !type
1131 ! Integrate analytic tail 1/(iw-mu)
1132  real(dp),intent(in) :: eig,fermie,temp
1133  real(dp) :: occup_fd
1134 !Local variables-------------------------------
1135 ! *********************************************************************
1136 
1137  if((eig-fermie) > zero) then
1138    occup_fd=exp(-(eig-fermie)/temp)/(one+exp(-(eig-fermie)/temp))
1139  else
1140    occup_fd=one/(one+exp((eig-fermie)/temp))
1141  endif
1142 
1143  end function occup_fd
1144 
1145 END MODULE m_energy

m_energy/print_energy [ Functions ]

[ Top ] [ m_energy ] [ Functions ]

NAME

 print_energy

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      m_energy

CHILDREN

      wrtout

SOURCE

251 subroutine print_energy(cryst_struc,energies_dmft,pawprtvol,pawtab,idmftloop)
252 
253 
254 !This section has been created automatically by the script Abilint (TD).
255 !Do not modify the following lines by hand.
256 #undef ABI_FUNC
257 #define ABI_FUNC 'print_energy'
258 !End of the abilint section
259 
260  implicit none
261 
262 !Arguments ------------------------------------
263 !type
264  type(crystal_t),intent(in) :: cryst_struc
265  type(energy_type),intent(in) :: energies_dmft
266  type(pawtab_type),intent(in)  :: pawtab(cryst_struc%ntypat)
267  integer, intent(in) :: pawprtvol,idmftloop
268 
269 !Local variables-------------------------------
270  integer :: iatom
271  character(len=1000) :: message
272 ! *********************************************************************
273  if(abs(pawprtvol)>=3) then
274    do iatom=1,cryst_struc%natom
275      if(pawtab(cryst_struc%typat(iatom))%lpawu/=-1) then
276        write(message,'(a,4x,a,i3,a,5x,f12.6)')  &
277 &       ch10,"For Correlated Atom",iatom,", E_hu =",energies_dmft%e_hu(iatom)
278        call wrtout(std_out,message,'COLL')
279        write(message,'(26x,a,1x,f12.6)')  &
280 &       ", E_hu_mig =",energies_dmft%e_hu_mig(iatom)
281        call wrtout(std_out,message,'COLL')
282        write(message,'(26x,a,f12.6)')  &
283 &       ", E_hu_qmc  =",energies_dmft%e_hu_qmc(iatom)
284        call wrtout(std_out,message,'COLL')
285        write(message,'(26x,a,f12.6)')  &
286 &       ", E_hu_ldau =",energies_dmft%e_hu_ldau(iatom)
287        call wrtout(std_out,message,'COLL')
288        write(message,'(26x,a,f12.6)')  &
289 &       ", E_dc =",energies_dmft%e_dc(iatom)
290        call wrtout(std_out,message,'COLL')
291      endif
292    enddo
293  endif
294  write(message,'(a,5x,2a,5x,a,9(a,5x,a,2x,f15.11),a,5x,a)') ch10 &
295 &      ,"-----------------------------------------------",ch10 &
296 &      ,"--- Energy in DMFT (in Ha)  ",ch10 &
297 &      ,"--- E_bandlda (1)  (Ha.) = ",energies_dmft%eband_lda,ch10 &
298 &      ,"--- E_banddmft(2)  (Ha.) = ",energies_dmft%eband_dmft,ch10 &
299 &      ,"--- E_hu      (3)  (Ha.) = ",energies_dmft%e_hu_tot,ch10 &
300 &      ,"--- E_hu_mig  (4)  (Ha.) = ",energies_dmft%e_hu_mig_tot,ch10 &
301 &      ,"--- E_hu_qmc  (4)  (Ha.) = ",energies_dmft%e_hu_qmc_tot,ch10 &
302 &      ,"--- E_hu_ldau (5)  (Ha.) = ",energies_dmft%e_hu_ldau_tot,ch10 &
303 &      ,"--- E_dc      (6)  (Ha.) = ",energies_dmft%e_dc_tot,ch10 &
304 &      ,"--- edmft=(    3-6)(Ha.) = ",energies_dmft%edmft,ch10 &
305 &      ,"---       (2-1+3-6)(Ha.) = ",energies_dmft%eband_dmft-energies_dmft%eband_lda+energies_dmft%edmft,ch10 &
306 &      ,"-----------------------------------------------"
307  call wrtout(std_out,message,'COLL')
308  if(idmftloop>=1) then
309    write(message,'(a,i3,1x,f15.11,a)') " (Edmft",idmftloop,energies_dmft%edmft,")"
310    call wrtout(ab_out,message,'COLL')
311  endif
312 end subroutine print_energy