TABLE OF CONTENTS
- ABINIT/m_energy
- m_energy/compute_band_energy
- m_energy/compute_dftu_energy
- m_energy/compute_energy
- m_energy/compute_migdal_energy
- m_energy/compute_noninterentropy
- m_energy/destroy_energy
- m_energy/energy_type
- m_energy/init_energy
- m_energy/occup_fd
- m_energy/print_energy
ABINIT/m_energy [ Modules ]
NAME
m_energy
FUNCTION
COPYRIGHT
Copyright (C) 2006-2024 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_energy 27 28 use defs_basis 29 use m_errors 30 use m_abicore 31 32 use m_pawtab, only : pawtab_type 33 use m_paw_correlations, only : pawuenergy 34 35 use m_green, only : green_type,icip_green,destroy_green,compa_occup_ks 36 use m_self, only : self_type,initialize_self,destroy_self,print_self,new_self,make_qmcshift_self 37 use m_paw_dmft, only : paw_dmft_type 38 39 use m_matlu, only : matlu_type,init_matlu,prod_matlu,diag_matlu,destroy_matlu,conjg_matlu,& 40 & ln_matlu,add_matlu,zero_matlu,shift_matlu,copy_matlu,trace_matlu,print_matlu 41 use m_oper, only : trace_oper 42 use m_crystal, only : crystal_t 43 44 implicit none 45 46 private 47 48 public :: init_energy 49 public :: compute_energy 50 public :: compute_migdal_energy 51 public :: compute_dftu_energy 52 public :: destroy_energy 53 public :: print_energy 54 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_dft= if present, compute free energy instead of total energy. fcalc_dft= 2 do a free energy calculation with the dmft fermie level = 2/3 subtract fermie to eigenvalues and free energy =1/3 fermie_dft level and free energy ecalc_dft= 1 subtract fermie to eigenvalues
OUTPUT
SIDE EFFECTS
energies_dmft <type(energy_type)> = DMFT energy structure data
SOURCE
455 subroutine compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_dft,fcalc_dft,ecalc_dmft) 456 457 !Arguments ------------------------------------ 458 !type 459 type(energy_type),intent(inout) :: energies_dmft 460 type(green_type),intent(in) :: green 461 type(paw_dmft_type), intent(inout) :: paw_dmft 462 character(len=4), intent(in) :: occ_type 463 integer, intent(in), optional :: ecalc_dft 464 integer, intent(in), optional :: fcalc_dft 465 integer, intent(in), optional :: ecalc_dmft 466 ! integer :: prtopt 467 468 !Local variables------------------------------- 469 integer :: ib,ikpt,isppol 470 real(dp) :: beta , fermie_used, totch,totch2,totch3 471 character(len=500) :: message 472 ! ********************************************************************* 473 write(message,'(2a)') ch10," == Compute Band Energy terms for DMFT " 474 call wrtout(std_out,message,'COLL') 475 beta=one/paw_dmft%temp 476 477 ! == Compute Band Energy 478 ! ----------------------------------------------------------------------- 479 energies_dmft%eband_dft=zero 480 energies_dmft%eband_dmft=zero 481 totch=zero 482 totch2=zero 483 totch3=zero 484 do isppol=1,paw_dmft%nsppol 485 do ikpt=1,paw_dmft%nkpt 486 do ib=1,paw_dmft%mbandc 487 if(present(fcalc_dft)) then 488 if (fcalc_dft==1.or.fcalc_dft==3) fermie_used=paw_dmft%fermie_dft 489 if (fcalc_dft==2.or.fcalc_dft==4) fermie_used=paw_dmft%fermie ! only for B3 terms 490 if((paw_dmft%eigen_dft(isppol,ikpt,ib)-fermie_used).ge.zero) then 491 energies_dmft%eband_dft=energies_dmft%eband_dft - & 492 & log ( one + exp ( - beta* (paw_dmft%eigen_dft(isppol,ikpt,ib)-fermie_used)))*paw_dmft%wtk(ikpt) 493 else 494 energies_dmft%eband_dft=energies_dmft%eband_dft & 495 & -(log ( one + exp ( beta* (paw_dmft%eigen_dft(isppol,ikpt,ib)-fermie_used))))*paw_dmft%wtk(ikpt) & 496 & + beta*paw_dmft%eigen_dft(isppol,ikpt,ib)*paw_dmft%wtk(ikpt) 497 if(fcalc_dft==3.or.fcalc_dft==2) then 498 ! subtract fermi level, (useful to directly count the number of electrons) 499 energies_dmft%eband_dft=energies_dmft%eband_dft & 500 & - beta*fermie_used*paw_dmft%wtk(ikpt) 501 totch=totch+paw_dmft%wtk(ikpt) 502 endif 503 endif 504 else ! usual calculation: total non interacting energy 505 fermie_used=paw_dmft%fermie_dft 506 ! write(std_out,*) "isppol,ikpt,ib",isppol,ikpt,ib 507 ! write(std_out,*) "paw_dmft%eigen_dft",paw_dmft%eigen_dft(isppol,ikpt,ib) 508 ! write(std_out,*) green%occup%ks(isppol,ikpt,ib,ib) 509 ! write(std_out,*) occup_fd(paw_dmft%eigen_dft(isppol,ikpt,ib),paw_dmft%fermie,paw_dmft%temp) 510 if(present(ecalc_dft)) then 511 if(ecalc_dft==1.or.ecalc_dft==3) fermie_used=paw_dmft%fermie_dft 512 if(ecalc_dft==2.or.ecalc_dft==4) fermie_used=paw_dmft%fermie ! only for B3 terms 513 if(ecalc_dft==3.or.ecalc_dft==2) then 514 energies_dmft%eband_dft=energies_dmft%eband_dft- & 515 & occup_fd(paw_dmft%eigen_dft(isppol,ikpt,ib),fermie_used,paw_dmft%temp)*& 516 & fermie_used*paw_dmft%wtk(ikpt) 517 totch2=totch2+paw_dmft%wtk(ikpt)*occup_fd(paw_dmft%eigen_dft(isppol,ikpt,ib),fermie_used,paw_dmft%temp) 518 endif 519 endif 520 energies_dmft%eband_dft=energies_dmft%eband_dft+ & 521 & occup_fd(paw_dmft%eigen_dft(isppol,ikpt,ib),fermie_used,paw_dmft%temp)*& 522 & paw_dmft%eigen_dft(isppol,ikpt,ib)*paw_dmft%wtk(ikpt) 523 endif 524 energies_dmft%eband_dmft=energies_dmft%eband_dmft+ & 525 & green%occup%ks(isppol,ikpt,ib,ib)*& 526 & paw_dmft%eigen_dft(isppol,ikpt,ib)*paw_dmft%wtk(ikpt) 527 totch3=totch3+paw_dmft%wtk(ikpt)*green%occup%ks(isppol,ikpt,ib,ib) 528 if(present(ecalc_dmft)) then 529 energies_dmft%eband_dmft=energies_dmft%eband_dmft- & 530 & green%occup%ks(isppol,ikpt,ib,ib)*& 531 & paw_dmft%fermie*paw_dmft%wtk(ikpt) 532 endif 533 enddo 534 enddo 535 enddo 536 if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) energies_dmft%eband_dft=two*energies_dmft%eband_dft 537 if(paw_dmft%nsppol==1.and.paw_dmft%nspinor==1) energies_dmft%eband_dmft=two*energies_dmft%eband_dmft 538 if(present(fcalc_dft)) then 539 energies_dmft%eband_dft=energies_dmft%eband_dft/beta 540 if(fcalc_dft==3.or.fcalc_dft==2) write(std_out,*) "compute_band_energy totch",totch 541 endif 542 if(present(ecalc_dft)) then 543 if(ecalc_dft==3.or.ecalc_dft==2) write(std_out,*) "compute_band_energy totch2",totch2 544 endif 545 ! write(std_out,*) "compute_band_energy totch3",totch3 546 547 if (occ_type==" lda") then 548 if(abs(energies_dmft%eband_dft-energies_dmft%eband_dmft)>tol5) then 549 write(message,'(5x,a,a,a,15x,a,f12.6,a,15x,a,5x,f12.5)') "Warning !:"& 550 & ,"Differences between band energy from DFT occupations",ch10& 551 & ,"and DFT green function is:",energies_dmft%eband_dft-energies_dmft%eband_dmft,ch10& 552 & ,"which is larger than",tol5 553 call wrtout(std_out,message,'COLL') 554 write(message,'(a)') & 555 & " Action: increase number of frequencies, or reduce the number of high energies_dmft bands" 556 call wrtout(std_out,message,'COLL') 557 else 558 write(message,'(a,a,a,10x,a,f12.6,a,10x,a,5x,f12.5)') " "& 559 & ,"Differences between band energy from DFT occupations",ch10& 560 & ,"and DFT green function is:",energies_dmft%eband_dft-energies_dmft%eband_dmft,ch10& 561 & ,"which is smaller than",tol5 562 call wrtout(std_out,message,'COLL') 563 endif 564 endif 565 566 567 end subroutine compute_band_energy
m_energy/compute_dftu_energy [ Functions ]
[ Top ] [ m_energy ] [ Functions ]
NAME
compute_dftu_energy
FUNCTION
Initialize noccmmp from green%occup and compute DFT+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
SOURCE
731 subroutine compute_dftu_energy(cryst_struc,energies_dmft,green,paw_dmft,pawtab,renorm) 732 733 !Arguments ------------------------------------ 734 !type 735 type(energy_type),intent(inout) :: energies_dmft 736 type(crystal_t),intent(in) :: cryst_struc 737 type(green_type),intent(in) :: green 738 type(paw_dmft_type), intent(in) :: paw_dmft 739 type(pawtab_type),target,intent(in) :: pawtab(cryst_struc%ntypat) 740 real(dp), optional, intent(in) :: renorm(:) 741 ! integer :: prtopt 742 743 !Local variables------------------------------- 744 integer :: iatom,idijeff,im,im1,ispinor,ispinor1,isppol,ldim,lpawu 745 integer :: nocc,nsploop,prt_pawuenergy 746 real(dp) :: upawu,jpawu 747 real(dp) :: edftumdcdc,edftumdc,e_ee,e_dc,e_dcdc,xe1,xe2 748 real(dp) :: edftumdc_for_s,edftumdcdc_for_s,e_ee_for_s,e_dc_for_s,e_dcdc_for_s 749 character(len=500) :: message 750 ! arrays 751 integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/)) 752 real(dp),allocatable :: noccmmp(:,:,:,:),nocctot(:) 753 type(pawtab_type),pointer :: pawtab_ 754 755 ! ********************************************************************* 756 757 ! - allocations 758 ! ----------------------------------------------------------------------- 759 760 nsploop=max(paw_dmft%nsppol,paw_dmft%nspinor**2) 761 nocc=nsploop 762 e_ee=zero 763 e_dc=zero 764 e_dc_for_s=zero 765 e_dcdc=zero 766 edftumdc=zero 767 edftumdcdc=zero 768 isppol=0 769 ispinor=0 770 ispinor1=0 771 772 ! - Loop and call to pawuenergy 773 ! ----------------------------------------------------------------------- 774 do iatom=1,cryst_struc%natom 775 pawtab_ => pawtab(cryst_struc%typat(iatom)) 776 lpawu=paw_dmft%lpawu(iatom) 777 if(lpawu.ne.-1) then 778 ldim=2*lpawu+1 779 780 ABI_MALLOC(noccmmp,(2,2*pawtab_%lpawu+1,2*pawtab_%lpawu+1,nocc)) 781 ABI_MALLOC(nocctot,(nocc)) 782 noccmmp(:,:,:,:)=zero ; nocctot(:)=zero 783 784 ! - Setup nocctot and noccmmp 785 ! ----------------------------------------------------------------------- 786 nocctot(:)=zero ! contains nmmp in the n m representation 787 ! Begin loop over spin/spinors to initialize noccmmp 788 do idijeff=1,nsploop 789 if(nsploop==2) then 790 isppol=spinor_idxs(1,idijeff) 791 ispinor=1 792 ispinor1=1 793 else if(nsploop==4) then 794 isppol=1 795 ispinor=spinor_idxs(1,idijeff) 796 ispinor1=spinor_idxs(2,idijeff) 797 else if(nsploop==1) then 798 isppol=1 799 ispinor=1 800 ispinor1=1 801 else 802 write(message,'(2a)') " BUG in m_energy: nsploop should be equal to 1, 2 or 4" 803 call wrtout(std_out,message,'COLL') 804 endif 805 ! Initialize noccmmp 806 do im1 = 1 , ldim 807 do im = 1 , ldim 808 noccmmp(1,im,im1,idijeff)=real(green%occup%matlu(iatom)%mat(im1,im,isppol,ispinor,ispinor1)) 809 noccmmp(2,im,im1,idijeff)=aimag(green%occup%matlu(iatom)%mat(im1,im,isppol,ispinor,ispinor1)) 810 ! noccmmp(1,im,im1,idijeff)=real(green%occup%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)) 811 ! noccmmp(2,im,im1,idijeff)=imag(green%occup%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)) 812 enddo 813 enddo 814 ! Compute nocctot 815 if(green%has_charge_matlu_solver/=2) then 816 do im1=1,ldim 817 if(nsploop==4) then 818 if(idijeff<=2) then 819 nocctot(1)=nocctot(1)+& 820 & noccmmp(1,im1,im1,idijeff) 821 endif 822 else if (nsploop==2.or.nsploop==1) then 823 nocctot(idijeff)=nocctot(idijeff)+& 824 & noccmmp(1,im1,im1,idijeff) 825 end if 826 enddo 827 else 828 if(nsploop==4) then 829 nocctot(1)=green%charge_matlu_solver(iatom,2) ! total nb of elec for nspinor=2 is (iatom,2) !! 830 nocctot(2)=zero 831 nocctot(3)=zero 832 nocctot(4)=zero 833 else if (nsploop==2) then 834 nocctot(1)=green%charge_matlu_solver(iatom,1) ! first spin 835 nocctot(2)=green%charge_matlu_solver(iatom,2) ! second one 836 else if (nsploop==1) then 837 nocctot(1)=green%charge_matlu_solver(iatom,1) 838 end if 839 endif 840 enddo 841 842 xe1=e_dc 843 xe2=e_ee 844 ! write(std_out,*)" nocctot(1)",nocctot(1),green%charge_matlu_solver(iatom,1) 845 edftumdc = zero 846 edftumdcdc = zero 847 if ( present(renorm) ) then 848 upawu = one 849 jpawu = renorm(iatom) 850 prt_pawuenergy=0 851 else 852 upawu = pawtab_%upawu 853 jpawu = pawtab_%jpawu 854 prt_pawuenergy=3 855 end if 856 857 call pawuenergy(iatom,edftumdc,edftumdcdc,noccmmp,nocctot,prt_pawuenergy,pawtab_,& 858 & dmft_dc=paw_dmft%dmft_dc,e_ee=e_ee,e_dc=e_dc,e_dcdc=e_dcdc,& 859 & u_dmft=upawu,j_dmft=jpawu) 860 861 if(paw_dmft%ientropy==1) then 862 call pawuenergy(iatom,edftumdc_for_s,edftumdcdc_for_s,noccmmp,nocctot,prt_pawuenergy,pawtab_,& 863 & dmft_dc=paw_dmft%dmft_dc,e_ee=e_ee_for_s,e_dc=e_dc_for_s,e_dcdc=e_dcdc_for_s,& 864 & u_dmft=paw_dmft%u_for_s/Ha_eV,j_dmft=paw_dmft%j_for_s/Ha_eV) 865 endif 866 867 868 energies_dmft%e_dc(iatom)=e_dc-xe1 ! probably wrong but not used 869 energies_dmft%e_hu_dftu(iatom)=e_ee-xe2 ! idem 870 871 ABI_FREE(noccmmp) 872 ABI_FREE(nocctot) 873 endif ! lpawu/=-1 874 enddo 875 876 ! - gather results 877 ! ----------------------------------------------------------------------- 878 energies_dmft%e_dc_tot=e_dc ! this is the onlu quantity used after. 879 energies_dmft%e_hu_dftu_tot=e_ee 880 if(paw_dmft%ientropy==1) then 881 write(message,'(a,3(f14.10,3x))') "For entropy calculation E_dc_tot, u_for_s, j_for,s", & 882 & e_dc_for_s,paw_dmft%u_for_s,paw_dmft%j_for_s 883 call wrtout(std_out,message,'COLL') 884 write(message,'(a,3(f14.10,3x))') "Reference calculation E_dc_tot, upawu , jpawu ", & 885 & e_dc,upawu*Ha_eV,jpawu*Ha_eV 886 call wrtout(std_out,message,'COLL') 887 endif 888 889 end subroutine compute_dftu_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
SOURCE
288 subroutine compute_energy(cryst_struc,energies_dmft,green,paw_dmft,pawprtvol,pawtab,self,occ_type,part) 289 290 !Arguments ------------------------------------ 291 !type 292 type(energy_type),intent(inout) :: energies_dmft 293 type(crystal_t),intent(in) :: cryst_struc 294 type(green_type),intent(in) :: green 295 type(paw_dmft_type), intent(inout) :: paw_dmft 296 type(pawtab_type),intent(in) :: pawtab(cryst_struc%ntypat) 297 type(self_type), intent(inout) :: self 298 integer, intent(in) :: pawprtvol 299 character(len=4), intent(in) :: occ_type 300 character(len=4), intent(in) :: part 301 ! integer :: prtopt 302 303 !Local variables------------------------------- 304 integer :: iatom,lpawu 305 integer :: natom,nspinor,nsppol 306 real(dp) :: beta 307 character(len=500) :: message 308 real(dp), allocatable :: e_hu_migdal(:) 309 real(dp) :: e_hu_migdal_tot 310 ! ********************************************************************* 311 if(part=='both') then 312 write(message,'(2a)') ch10," == Compute DFT+DMFT energy terms " 313 call wrtout(std_out,message,'COLL') 314 else if(part=='band') then 315 write(message,'(2a)') ch10," == Compute DFT+DMFT energy terms : Band energy terms" 316 call wrtout(std_out,message,'COLL') 317 else if(part=='corr') then 318 write(message,'(2a)') ch10," == Compute DFT+DMFT energy terms : Correlation energy terms only" 319 call wrtout(std_out,message,'COLL') 320 else if(part=='none') then 321 endif 322 323 ! Only imaginary frequencies here 324 if(green%w_type=="real".or.self%w_type=="real") then 325 message = 'compute_energy not implemented for real frequency' 326 ABI_BUG(message) 327 endif 328 natom=cryst_struc%natom 329 nsppol=paw_dmft%nsppol 330 nspinor=paw_dmft%nspinor 331 beta=one/paw_dmft%temp 332 333 if(part=='band'.or.part=='both') then 334 ! == Compute Band Energy Alternative version: two steps 335 ! == Compute Tr[ln G^{-1}] and -Tr[(Self-hdc)G_dmft] 336 ! ----------------------------------------------------------------------- 337 if(part=='band') then ! ie if thdyn="fcalc" in m_dmft.F90 338 ! call compute_B3(cryst_struc,energies_dmft,eband2,green,mpi_enreg,paw_dmft,2,pawang,self,occ_type,0) 339 ! write(message,'(2a,f10.6)') ch10,"Compute Band energy test KS ",eband2 340 ! call wrtout(std_out,message,'COLL') 341 ! call wrtout(ab_out,message,'COLL') 342 ! 343 ! call compute_B3(cryst_struc,energies_dmft,eband2,green,mpi_enreg,paw_dmft,2,pawang,self,occ_type,1) 344 ! write(message,'(2a,f10.6)') ch10,"Compute Band energy test Self statique",eband2 345 ! call wrtout(std_out,message,'COLL') 346 ! call wrtout(ab_out,message,'COLL') 347 ! 348 !! == Compute Band Energy (classical) 349 !! ----------------------------------------------------------------------- 350 ! call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,fcalc_dft=3) 351 ! write(std_out,*) paw_dmft%fermie_dft,paw_dmft%fermie 352 ! write(message,'(2a,f10.6)') ch10,"Compute Band energy ref free lda -ef ",energies_dmft%eband_dft 353 ! call wrtout(std_out,message,'COLL') 354 call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_dft=1) 355 ! write(message,'(2a,f10.6)') ch10,"Compute Band energy ref -ef ",energies_dmft%eband_dmft 356 ! call wrtout(std_out,message,'COLL') 357 !! call wrtout(ab_out,message,'COLL') 358 ! write(message,'(2a,f10.6)') ch10,"Compute Band energy ref lda ",energies_dmft%eband_dft 359 ! call wrtout(std_out,message,'COLL') 360 !! if(occ_type=="nlda") eband2=energies_dmft%eband_dmft 361 else 362 call compute_band_energy(energies_dmft,green,paw_dmft,occ_type,ecalc_dft=1) 363 endif 364 365 endif 366 if(part=='corr'.or.part=='both') then 367 368 ! == Compute Correlation energy from Migdal formula 369 ! ----------------------------------------------------------------------- 370 ABI_MALLOC(e_hu_migdal,(cryst_struc%natom)) 371 e_hu_migdal(:) = zero 372 call compute_migdal_energy(cryst_struc,e_hu_migdal,e_hu_migdal_tot,green,paw_dmft,pawprtvol,self) 373 energies_dmft%e_hu_mig(:)= e_hu_migdal(:) 374 energies_dmft%e_hu_mig_tot = e_hu_migdal_tot 375 ! write(std_out,*) "MIGDAL",e_hu_migdal_tot,e_hu_migdal 376 ABI_FREE(e_hu_migdal) 377 378 ! == Compute Correlation energy from QMC correlations. 379 ! ----------------------------------------------------------------------- 380 energies_dmft%e_hu_qmc_tot = zero 381 do iatom=1,natom 382 lpawu=paw_dmft%lpawu(iatom) 383 if(lpawu/=-1) then 384 if(paw_dmft%dmft_solv==4.or.paw_dmft%dmft_solv==5.or.paw_dmft%dmft_solv==8) then 385 energies_dmft%e_hu_qmc(iatom) = green%ecorr_qmc(iatom) 386 energies_dmft%e_hu_qmc_tot = energies_dmft%e_hu_qmc_tot + green%ecorr_qmc(iatom) 387 endif 388 endif ! lpawu 389 enddo ! iatom 390 391 ! == Compute DFT+U interaction energy 392 ! ----------------------------------------------------------------------- 393 call compute_dftu_energy(cryst_struc,energies_dmft,green,paw_dmft,pawtab) 394 if(abs(paw_dmft%dmft_solv)<=1) then 395 energies_dmft%e_hu= energies_dmft%e_hu_dftu 396 energies_dmft%e_hu_tot= energies_dmft%e_hu_dftu_tot 397 if((abs(energies_dmft%e_hu_tot-energies_dmft%e_hu_mig_tot).ge.0.000001).and.(occ_type/=" lda")) then 398 write(message,'(2a,2e18.8,a)') ch10,' BUG: Migdal energy and DFT+U energy do not coincide',& 399 & energies_dmft%e_hu_tot,energies_dmft%e_hu_mig_tot,occ_type 400 ABI_ERROR(message) 401 endif 402 else if(paw_dmft%dmft_solv==2.or.paw_dmft%dmft_solv==6.or.paw_dmft%dmft_solv==7.or.paw_dmft%dmft_solv==9) then 403 energies_dmft%e_hu= energies_dmft%e_hu_mig 404 energies_dmft%e_hu_tot= energies_dmft%e_hu_mig_tot 405 energies_dmft%e_hu_qmc_tot = energies_dmft%e_hu_tot 406 else if(paw_dmft%dmft_solv==4.or.paw_dmft%dmft_solv==5.or.paw_dmft%dmft_solv==8) then 407 if(paw_dmft%dmft_solv==8) then 408 write(message,'(2a)') ch10,"Warning, energy is recently computed, not checked" 409 call wrtout(std_out,message,'COLL') 410 endif 411 energies_dmft%e_hu= energies_dmft%e_hu_qmc 412 energies_dmft%e_hu_tot= energies_dmft%e_hu_qmc_tot 413 endif 414 ! energies_dmft%edmft=energies_dmft%e_hu_mig_tot-energies_dmft%e_dc_tot 415 energies_dmft%edmft=energies_dmft%e_hu_tot-energies_dmft%e_dc_tot 416 417 418 endif ! part 419 420 ! if(part='corr'.or.part='both') then 421 if(part/='none') then 422 call print_energy(cryst_struc,energies_dmft,pawprtvol,pawtab,paw_dmft%idmftloop) 423 endif 424 ! write(message,'(2a)') ch10," == The DFT+U self-energy is == " 425 ! call wrtout(std_out,message,'COLL') 426 ! call print_oper(self%oper(1),5,paw_dmft,2) 427 ! a voir: energies_dmft%e_hu_tot = energies_dmft%e_hu_dftu_tot 428 429 end subroutine compute_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.
SOURCE
591 subroutine compute_migdal_energy(cryst_struc,e_hu_migdal,e_hu_migdal_tot,green,paw_dmft,pawprtvol,self) 592 593 #ifdef FC_INTEL 594 !DEC$ NOOPTIMIZE 595 #endif 596 597 !Arguments ------------------------------------ 598 !type 599 type(crystal_t),intent(in) :: cryst_struc 600 type(green_type),intent(in) :: green 601 type(paw_dmft_type), intent(in) :: paw_dmft 602 real(dp),intent(out) :: e_hu_migdal_tot 603 real(dp),intent(out) :: e_hu_migdal(paw_dmft%natom) 604 type(self_type), intent(in) :: self 605 integer, intent(in) :: pawprtvol 606 ! integer :: prtopt 607 608 !Local variables------------------------------- 609 integer :: iatom,ifreq,im,im1,ispinor,ispinor1,isppol,lpawu 610 integer :: natom,ndim,nspinor,nsppol,nwlo 611 real(dp) :: beta 612 complex(dpc) :: xmig_1,xmig_2,xmig_3,se,shift 613 character(len=500) :: message 614 ! ********************************************************************* 615 616 ! Only imaginary frequencies here 617 if(green%w_type=="real".or.self%w_type=="real") then 618 message = 'compute_migdal_energy not implemented for real frequency' 619 ABI_BUG(message) 620 endif 621 622 ! == Compute Correlation energy from Migdal formula 623 ! ----------------------------------------------------------------------- 624 natom=cryst_struc%natom 625 nsppol=paw_dmft%nsppol 626 nspinor=paw_dmft%nspinor 627 beta=one/paw_dmft%temp 628 nwlo=green%nw 629 if (green%nw/=self%nw) then 630 message = 'self and green do not contains the same number of frequencies' 631 ABI_BUG(message) 632 endif 633 ! write(std_out,*) "beta",beta 634 635 xmig_1=zero 636 xmig_2=zero 637 xmig_3=zero 638 639 e_hu_migdal_tot = zero 640 do iatom=1,natom 641 shift=czero 642 if(paw_dmft%dmft_solv==4) shift=self%qmc_shift(iatom)+self%qmc_xmu(iatom) 643 ! write(std_out,*) "shiftttt",shift 644 lpawu=paw_dmft%lpawu(iatom) 645 if(lpawu/=-1) then 646 xmig_1=czero 647 xmig_2=czero 648 xmig_3=czero 649 ndim=2*lpawu+1 650 do isppol=1,nsppol 651 do ispinor = 1 , nspinor 652 do ispinor1 = 1, nspinor 653 do im=1,ndim 654 do im1=1,ndim 655 do ifreq=1,nwlo 656 ! write(std_out,*) ifreq,xmig_1,imag(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)),& 657 !& green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor ) 658 xmig_1=xmig_1 + j_dpc/beta* & 659 & aimag(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1))* & 660 & green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor )* & 661 & paw_dmft%wgt_wlo(ifreq) 662 ! if(ispinor==ispinor1.and.im==im1) then 663 se=(self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)- & 664 & self%oper (nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)) 665 ! else 666 ! se=self%oper (ifreq)%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1) 667 ! endif 668 xmig_2=xmig_2 + one/beta*real(se)* & 669 & green%oper(ifreq)%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor )* & 670 & paw_dmft%wgt_wlo(ifreq) 671 ! if(ispinor==ispinor1.nd.im==im1.and.ifreq==1) then 672 if(ifreq==1) then 673 xmig_3=xmig_3 + & 674 & real(self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1)+shift)* & 675 & green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two 676 ! write(std_out,*) "xmig_3",xmig_3 677 ! write(std_out,*) "self",self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1) 678 ! write(std_out,*) "shift",shift 679 ! write(std_out,*) "occup", green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two 680 endif 681 enddo 682 ! if(ispinor==ispinor1.and.im==im1) then 683 ! xmig_3=xmig_3 + & 684 !& real(self%oper(nwlo )%matlu(iatom)%mat(im ,im1,isppol,ispinor ,ispinor1))* & 685 !!& green%occup%matlu(iatom)%mat(im1,im ,isppol,ispinor1,ispinor)/two 686 ! endif 687 enddo 688 enddo 689 enddo 690 enddo 691 enddo 692 if(nsppol==1.and.nspinor==1) then 693 e_hu_migdal(iatom)=two*real(xmig_1+xmig_2+xmig_3) 694 else 695 e_hu_migdal(iatom)=real(xmig_1+xmig_2+xmig_3) 696 endif 697 e_hu_migdal_tot = e_hu_migdal_tot + e_hu_migdal(iatom) 698 if(abs(pawprtvol)>=3) then 699 write(message,'(2a,3(a,5x,a,2f12.6))')ch10,& 700 & " Interaction energy: Decomposition of Migdal energy",ch10,& 701 & "xmig_1=",xmig_1,ch10,& 702 & "xmig_2=",xmig_2,ch10,& 703 & "xmig_3=",xmig_3 704 call wrtout(std_out,message,'COLL') 705 endif 706 endif ! lpawu 707 enddo 708 709 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
SOURCE
909 subroutine compute_noninterentropy(cryst_struc,green,paw_dmft) 910 911 !Arguments ------------------------------------ 912 !type 913 type(crystal_t),intent(in) :: cryst_struc 914 type(green_type),intent(in) :: green 915 type(paw_dmft_type), intent(inout) :: paw_dmft 916 917 !Local variables------------------------------- 918 integer :: ib,ikpt,isppol,natom,nspinor,nsppol 919 real(dp) :: beta,eig,fermi,s_1,s_2,occ1,occ2,f_1,e_1,f_1a,s_1a,e_2 920 character(len=800) :: message 921 ! ********************************************************************* 922 write(message,'(2a,i6)') ch10," == Compute T*Entropy for fermi level and DFT-KS eigenvalues " 923 call wrtout(std_out,message,'COLL') 924 925 natom=cryst_struc%natom 926 nsppol=paw_dmft%nsppol 927 nspinor=paw_dmft%nspinor 928 beta=one/paw_dmft%temp 929 s_1=zero 930 s_1a=zero 931 f_1=zero 932 f_1a=zero 933 e_1=zero 934 e_2=zero 935 s_2=zero 936 do isppol=1,paw_dmft%nsppol 937 do ikpt=1,paw_dmft%nkpt 938 do ib=1,paw_dmft%mbandc 939 eig=paw_dmft%eigen_dft(isppol,ikpt,ib) 940 fermi=paw_dmft%fermie_dft 941 fermi=paw_dmft%fermie 942 occ1=occup_fd(eig,fermi,paw_dmft%temp) 943 occ2=green%occup%ks(isppol,ikpt,ib,ib) 944 ! write(std_out,*) occ1,occ2 945 946 ! entropy from Fermi Dirac 947 if((occ1.ge.tol9).and.((one-occ1).ge.tol9)) then 948 s_1=s_1+(occ1*log(occ1)+(one-occ1)*log(one-occ1))*paw_dmft%wtk(ikpt) 949 ! write(std_out,*) occ1,one-occ1,"p1",(occ1*log(occ1)+(one-occ1)*log(one-occ1))*paw_dmft%wtk(ikpt) 950 endif 951 952 ! Free energy from Fermi Dirac 953 if((eig-fermi).ge.zero) then ! occ1 -> 0 ; 1-occ1 -> 1 954 f_1=f_1-paw_dmft%wtk(ikpt)/beta*log(one+exp(-beta*(eig-fermi))) 955 f_1a=f_1a+paw_dmft%wtk(ikpt)/beta*log(one-occ1) 956 s_1a=s_1a+((one-occ1)*log(one-occ1)+occ1*(-beta*(eig-fermi)+log(one-occ1)))*paw_dmft%wtk(ikpt) 957 else ! occ1 -> 1 , 1-occ1 -> 0 958 f_1=f_1-paw_dmft%wtk(ikpt)/beta*(log(one+exp(beta*(eig-fermi)))-beta*(eig-fermi)) 959 f_1a=f_1a-paw_dmft%wtk(ikpt)/beta*(-log(occ1)-beta*(eig-fermi)) 960 s_1a=s_1a+(occ1*log(occ1)+(one-occ1)*(beta*(eig-fermi)+log(occ1)))*paw_dmft%wtk(ikpt) 961 endif 962 963 ! Internal energy from Fermi Dirac 964 e_1=e_1+(eig-fermi)*paw_dmft%wtk(ikpt)*occ1 965 e_2=e_2+(eig-fermi)*paw_dmft%wtk(ikpt)*occ2 966 967 ! entropy from green function occupations. 968 if((occ2.ge.tol9).and.((one-occ2).ge.tol9)) then 969 ! write(std_out,*) occ2,one-occ2,"p2",(occ2*log(occ2)+(one-occ2)*log(one-occ2))*paw_dmft%wtk(ikpt) 970 s_2=s_2+(occ2*log(occ2)+(one-occ2)*log(one-occ2))*paw_dmft%wtk(ikpt) 971 endif 972 enddo 973 enddo 974 enddo 975 s_1=-s_1*paw_dmft%temp 976 s_1a=-s_1a*paw_dmft%temp 977 s_2=-s_2*paw_dmft%temp 978 979 write(message,'(8(2a,e20.9))') & 980 & ch10," T*Entropy from Fermi Dirac occupations ", s_1,& 981 & ch10," T*Entropy from Fermi Dirac occupations 2 ", s_1a,& 982 & ch10," T*Entropy from Green function occupations ", s_2,& 983 & ch10," Free energy F ", f_1,& 984 & ch10," Free energy Fa ", f_1a,& 985 & ch10," internal energy U ", e_1,& 986 & ch10," internal energy U from Gr Func Occ ", e_2,& 987 & ch10," U-F ", e_1-f_1 988 call wrtout(std_out,message,'COLL') 989 990 991 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
SOURCE
167 subroutine destroy_energy(energies_dmft,paw_dmft) 168 169 !Arguments ------------------------------------ 170 !scalars 171 type(energy_type),intent(inout) :: energies_dmft 172 type(paw_dmft_type), intent(inout) :: paw_dmft 173 !Local variables------------------------------- 174 ! ********************************************************************* 175 paw_dmft%edmft=energies_dmft%edmft 176 if ( allocated(energies_dmft%e_dc) ) then 177 ABI_FREE(energies_dmft%e_dc) 178 end if 179 if ( allocated(energies_dmft%e_hu) ) then 180 ABI_FREE(energies_dmft%e_hu) 181 end if 182 if ( allocated(energies_dmft%e_hu_dftu) ) then 183 ABI_FREE(energies_dmft%e_hu_dftu) 184 end if 185 if ( allocated(energies_dmft%e_hu_mig) ) then 186 ABI_FREE(energies_dmft%e_hu_mig) 187 end if 188 if ( allocated(energies_dmft%e_hu_qmc) ) then 189 ABI_FREE(energies_dmft%e_hu_qmc) 190 end if 191 192 193 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
67 type, public :: energy_type ! for each typat 68 69 real(dp) :: eband_dft 70 71 real(dp) :: eband_dmft 72 73 real(dp) :: e_dc_tot 74 75 real(dp) :: e_hu_tot 76 77 real(dp) :: e_hu_dftu_tot 78 79 real(dp) :: e_hu_mig_tot 80 81 real(dp) :: e_hu_qmc_tot 82 83 real(dp) :: edmft 84 85 real(dp) :: natom 86 87 real(dp), allocatable :: e_dc(:) 88 89 real(dp), allocatable :: e_hu(:) 90 91 real(dp), allocatable :: e_hu_dftu(:) 92 93 real(dp), allocatable :: e_hu_mig(:) 94 95 real(dp), allocatable :: e_hu_qmc(:) 96 97 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
SOURCE
121 subroutine init_energy(cryst_struc,energies_dmft) 122 123 !Arguments ------------------------------------ 124 !type 125 type(crystal_t),intent(in) :: cryst_struc 126 type(energy_type), intent(inout) :: energies_dmft 127 !Local variables ------------------------------------ 128 !************************************************************************ 129 130 ABI_MALLOC(energies_dmft%e_dc,(cryst_struc%natom)) 131 ABI_MALLOC(energies_dmft%e_hu,(cryst_struc%natom)) 132 ABI_MALLOC(energies_dmft%e_hu_dftu,(cryst_struc%natom)) 133 ABI_MALLOC(energies_dmft%e_hu_mig,(cryst_struc%natom)) 134 ABI_MALLOC(energies_dmft%e_hu_qmc,(cryst_struc%natom)) 135 energies_dmft%e_dc=zero 136 energies_dmft%e_hu=zero 137 energies_dmft%e_hu_dftu=zero 138 energies_dmft%e_hu_mig=zero 139 energies_dmft%e_hu_qmc=zero 140 energies_dmft%eband_dft=zero 141 energies_dmft%eband_dmft=zero 142 energies_dmft%e_dc_tot=zero 143 energies_dmft%e_hu_tot=zero 144 energies_dmft%e_hu_dftu_tot=zero 145 energies_dmft%e_hu_mig_tot=zero 146 energies_dmft%e_hu_qmc_tot=zero 147 energies_dmft%edmft=zero 148 energies_dmft%natom=cryst_struc%natom 149 150 end subroutine init_energy
m_energy/occup_fd [ Functions ]
[ Top ] [ m_energy ] [ Functions ]
NAME
occup_fd
FUNCTION
INPUTS
OUTPUT
SOURCE
1007 function occup_fd(eig,fermie,temp) 1008 1009 !Arguments ------------------------------------ 1010 !type 1011 ! Integrate analytic tail 1/(iw-mu) 1012 real(dp),intent(in) :: eig,fermie,temp 1013 real(dp) :: occup_fd 1014 !Local variables------------------------------- 1015 ! ********************************************************************* 1016 1017 if((eig-fermie) > zero) then 1018 occup_fd=exp(-(eig-fermie)/temp)/(one+exp(-(eig-fermie)/temp)) 1019 else 1020 occup_fd=one/(one+exp((eig-fermie)/temp)) 1021 endif 1022 1023 end function occup_fd 1024 1025 END MODULE m_energy
m_energy/print_energy [ Functions ]
[ Top ] [ m_energy ] [ Functions ]
NAME
print_energy
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
SOURCE
210 subroutine print_energy(cryst_struc,energies_dmft,pawprtvol,pawtab,idmftloop) 211 212 !Arguments ------------------------------------ 213 !type 214 type(crystal_t),intent(in) :: cryst_struc 215 type(energy_type),intent(in) :: energies_dmft 216 type(pawtab_type),intent(in) :: pawtab(cryst_struc%ntypat) 217 integer, intent(in) :: pawprtvol,idmftloop 218 219 !Local variables------------------------------- 220 integer :: iatom 221 character(len=1000) :: message 222 ! ********************************************************************* 223 if(abs(pawprtvol)>=3) then 224 do iatom=1,cryst_struc%natom 225 if(pawtab(cryst_struc%typat(iatom))%lpawu/=-1) then 226 write(message,'(a,4x,a,i3,a,5x,f12.6)') & 227 & ch10,"For Correlated Atom",iatom,", E_hu =",energies_dmft%e_hu(iatom) 228 call wrtout(std_out,message,'COLL') 229 write(message,'(26x,a,1x,f12.6)') & 230 & ", E_hu_mig =",energies_dmft%e_hu_mig(iatom) 231 call wrtout(std_out,message,'COLL') 232 write(message,'(26x,a,f12.6)') & 233 & ", E_hu_qmc =",energies_dmft%e_hu_qmc(iatom) 234 call wrtout(std_out,message,'COLL') 235 write(message,'(26x,a,f12.6)') & 236 & ", E_hu_dftu =",energies_dmft%e_hu_dftu(iatom) 237 call wrtout(std_out,message,'COLL') 238 write(message,'(26x,a,f12.6)') & 239 & ", E_dc =",energies_dmft%e_dc(iatom) 240 call wrtout(std_out,message,'COLL') 241 endif 242 enddo 243 endif 244 write(message,'(a,5x,2a,5x,a,9(a,5x,a,2x,f15.11),a,5x,a)') ch10 & 245 & ,"-----------------------------------------------",ch10 & 246 & ,"--- Energy in DMFT (in Ha) ",ch10 & 247 & ,"--- E_bandlda (1) (Ha.) = ",energies_dmft%eband_dft,ch10 & 248 & ,"--- E_banddmft(2) (Ha.) = ",energies_dmft%eband_dmft,ch10 & 249 & ,"--- E_hu (3) (Ha.) = ",energies_dmft%e_hu_tot,ch10 & 250 & ,"--- E_hu_mig (4) (Ha.) = ",energies_dmft%e_hu_mig_tot,ch10 & 251 & ,"--- E_hu_qmc (4) (Ha.) = ",energies_dmft%e_hu_qmc_tot,ch10 & 252 & ,"--- E_hu_dftu (5) (Ha.) = ",energies_dmft%e_hu_dftu_tot,ch10 & 253 & ,"--- E_dc (6) (Ha.) = ",energies_dmft%e_dc_tot,ch10 & 254 & ,"--- edmft=( 3-6)(Ha.) = ",energies_dmft%edmft,ch10 & 255 & ,"--- (2-1+3-6)(Ha.) = ",energies_dmft%eband_dmft-energies_dmft%eband_dft+energies_dmft%edmft,ch10 & 256 & ,"-----------------------------------------------" 257 call wrtout(std_out,message,'COLL') 258 if(idmftloop>=1) then 259 write(message,'(a,i3,1x,f15.11,a)') " (Edmft",idmftloop,energies_dmft%edmft,")" 260 call wrtout(ab_out,message,'COLL') 261 endif 262 end subroutine print_energy