TABLE OF CONTENTS


ABINIT/m_hubbard_one [ Modules ]

[ Top ] [ Modules ]

NAME

  m_hubbard_one

FUNCTION

 Solve Anderson model with the density/density Hubbard one approximation

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

25 #if defined HAVE_CONFIG_H
26 #include "config.h"
27 #endif
28 
29 
30 #include "abi_common.h"
31 
32 MODULE m_hubbard_one
33 
34 
35  use defs_basis
36 
37  implicit none
38 
39  private
40 
41  public :: hubbard_one

m_hubbard_one/combin [ Functions ]

[ Top ] [ m_hubbard_one ] [ Functions ]

NAME

 combin

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

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

811  recursive subroutine combin(ielec,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup)
812 
813  use defs_basis
814 
815 !This section has been created automatically by the script Abilint (TD).
816 !Do not modify the following lines by hand.
817 #undef ABI_FUNC
818 #define ABI_FUNC 'combin'
819 !End of the abilint section
820 
821  implicit none
822 
823 !Arguments ------------------------------------
824 !scalars
825 ! type(pawang_type), intent(in) :: pawang
826  integer, intent(in) :: ielec,nelec,nlevels
827  integer, intent(inout) :: nconfig,nconfig_nelec(0:nlevels)
828  integer, intent(inout) :: occup(0:nlevels,nlevels)
829 ! type  :: level2_type
830 !  integer, pointer :: repart(:,:)
831 ! end type
832  type(level2_type), intent(inout) :: occ_level(0:nlevels)
833 ! integer, intent(in) :: prtopt
834 
835 !Local variables ------------------------------
836 ! scalars
837  integer :: max_ielec,pos,min_ielec,jelec,prtopt
838  character(len=500) :: message
839 ! arrays
840 !************************************************************************
841    prtopt=1
842    max_ielec=nlevels-nelec+ielec
843 !  write(std_out,*) "call to combin ielec,nelec,nlevels",ielec,nelec,nlevels
844    select case (ielec)
845    case (1)
846      min_ielec=1
847    case default
848      min_ielec=occup(nelec,ielec-1)+1
849    end select
850 !  write(std_out,*) "For ielec", ielec, "min_ielec,max_ielec",min_ielec,max_ielec
851    do pos = min_ielec, max_ielec
852      if(ielec==nelec) then
853        occup(nelec,ielec)=pos
854        nconfig=nconfig+1
855        nconfig_nelec(nelec)=nconfig_nelec(nelec)+1
856 !      write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,1)
857 !      write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,2)
858        do jelec=1,nelec
859 !        write(std_out,*) "nconfig",nconfig_nelec(nelec),nelec
860 !        write(std_out,*) "occup",occup(nelec,jelec)
861          occ_level(nelec)%repart(nconfig_nelec(nelec),jelec)=occup(nelec,jelec)
862        end do
863 !      write(std_out,*) "For ielec", ielec, "case nelec"
864        if(prtopt>=3) then
865          write(message,'(a,i3,a,30i5)') "For ielec",ielec," Occupf are", (occup(nelec,jelec),jelec=1,nelec)
866          call wrtout(std_out,message,'COLL')
867        end if
868      else
869        occup(nelec,ielec)=pos
870 !      write(std_out,*) "For ielec", ielec, "case 1 and default"
871        call combin(ielec+1,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup)
872        if(prtopt>=3) then
873          write(message,'(a,i3,a,30i5)') "For ielec",ielec," Occup are", (occup(nelec,jelec),jelec=1,nelec)
874          call wrtout(std_out,message,'COLL')
875        end if
876      end if
877    end do
878 
879  end subroutine combin

m_hubbard_one/green_atomic_hubbard [ Functions ]

[ Top ] [ m_hubbard_one ] [ Functions ]

NAME

 green_atomic_hubbard

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

  cryst_struc
  istep    =  step of iteration for LDA.
  lda_occup
  mpi_enreg=informations about MPI parallelization
  paw_dmft =  data for self-consistent LDA+DMFT calculations.
  pawang <type(pawang)>=paw angular mesh and related data

OUTPUT

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

NOTES

PARENTS

      hubbard_one

CHILDREN

      combin,destroy_green,init_green,wrtout

SOURCE

402 subroutine green_atomic_hubbard(cryst_struc,green_hubbard,hu,level_diag,paw_dmft,udens_atoms)
403 
404 
405  use defs_basis
406  use m_errors
407  use m_abicore
408  use m_crystal, only : crystal_t
409  use m_special_funcs,  only : factorial, permutations
410  use m_green, only : green_type,init_green,destroy_green
411  use m_hu, only : hu_type
412  use m_paw_dmft, only : paw_dmft_type
413 
414 !This section has been created automatically by the script Abilint (TD).
415 !Do not modify the following lines by hand.
416 #undef ABI_FUNC
417 #define ABI_FUNC 'green_atomic_hubbard'
418 !End of the abilint section
419 
420  implicit none
421 
422 !Arguments ------------------------------------
423 !scalars
424 ! type(pawang_type), intent(in) :: pawang
425  type(crystal_t),intent(in) :: cryst_struc
426  type(green_type), intent(inout) :: green_hubbard
427  type(paw_dmft_type), intent(in)  :: paw_dmft
428  type(matlu_type),intent(in) :: level_diag(cryst_struc%natom)
429  type(hu_type), intent(inout) :: hu(cryst_struc%ntypat)
430  type(coeff2_type),intent(in) :: udens_atoms(cryst_struc%natom)
431 
432 !Local variables ------------------------------
433 ! scalars
434  integer :: cnk,iacc,iatom,iconfig,ielec,ifreq,ilevel,im,im1,isppol,ispinor,itrans,jconfig,jelec
435  integer :: lpawu,m_temp,nconfig,nelec,nlevels,nspinor,nsppol,occupied_level,sum_test
436  integer, allocatable :: occup(:,:),nconfig_nelec(:)
437  character(len=500) :: message
438 ! arrays
439  type(green_type) :: green_hubbard_realw
440  type(level2_type), allocatable :: occ_level(:)
441  type(level1_type), allocatable :: e_nelec(:)
442  complex(dpc), allocatable :: green_temp(:,:)
443  complex(dpc), allocatable :: green_temp_realw(:,:)
444  complex(dpc) :: Z_part
445  real(dp), allocatable :: maxener(:),minener(:)
446  real(dp), allocatable :: elevels(:)
447  real(dp) :: emax,emin,eshift,prtopt, Ej_np1, Ei_n,beta,maxarg_exp,tmp
448 !************************************************************************
449    maxarg_exp=300
450 
451 !  hu is not used anymore.
452    if(hu(1)%lpawu==0) then
453    end if
454 !  ======================================
455 !  General loop over atoms
456 !  ======================================
457    nsppol=paw_dmft%nsppol
458    nspinor=paw_dmft%nspinor
459    prtopt=1
460    beta=one/paw_dmft%temp
461    call init_green(green_hubbard_realw,paw_dmft,opt_oper_ksloc=2,wtype=green_hubbard%w_type) ! initialize only matlu
462 
463    do iatom=1,cryst_struc%natom
464      lpawu=paw_dmft%lpawu(iatom)
465      if(lpawu/=-1) then
466        nlevels=nsppol*nspinor*(2*lpawu+1)
467        if(nsppol==1.and.nspinor==1) nlevels=2*nspinor*(2*lpawu+1)
468 
469 !      ===================================
470 !      Allocations
471 !      ===================================
472        ABI_DATATYPE_ALLOCATE(occ_level,(0:nlevels))
473        ABI_ALLOCATE(maxener,(0:nlevels))
474        ABI_ALLOCATE(minener,(0:nlevels))
475        ABI_ALLOCATE(elevels,(nlevels))
476        ABI_DATATYPE_ALLOCATE(e_nelec,(0:nlevels))
477        do nelec=0,nlevels ! number of electrons
478          cnk=nint(permutations(nlevels,nelec)/factorial(nelec))
479          ABI_ALLOCATE(occ_level(nelec)%repart      ,(cnk,nelec))
480          ABI_ALLOCATE(occ_level(nelec)%ocp         ,(cnk,nlevels))
481          ABI_ALLOCATE(occ_level(nelec)%transition  ,(cnk,nlevels-nelec))
482          ABI_ALLOCATE(occ_level(nelec)%transition_m,(cnk,nlevels))
483          ABI_ALLOCATE(e_nelec  (nelec)%config      ,(cnk))
484          e_nelec(nelec)%config(:)=zero
485 !        write(std_out,*) "permutations",nint(permutations(nlevels,nelec)/factorial(nelec))
486 !        write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,1)
487 !        write(std_out,*) "size",size(occ_level),size(occ_level(nelec)%repart,2)
488 !        for a given nb of electrons nelec, gives for a given repartition
489 !        of electron, the position of the ielec electron inside atomic
490 !        levels
491 !        levels
492        end do
493        ABI_ALLOCATE(occup,(0:nlevels,nlevels))
494        ABI_ALLOCATE(nconfig_nelec,(0:nlevels))
495 
496 !      ===================================
497 !      Initialization
498 !      ===================================
499        nconfig_nelec=0
500        nconfig=1
501        occup=0
502        nconfig_nelec(0)=1
503        occup(0,:)=0
504        iacc=0
505        elevels=zero
506        do isppol=1,nsppol
507          do ispinor=1,nspinor
508            do im1=1,(2*lpawu+1)
509              iacc=iacc+1
510              elevels(iacc)=level_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor)
511              if(abs(aimag(level_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor)))>tol8) then
512                message = " Hubbard I: levels are imaginary"
513                write(std_out,*) level_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor)
514                MSG_BUG(message)
515              end if
516              if(nsppol==1.and.nspinor==1) then
517                elevels(nspinor*(2*lpawu+1)+iacc)=level_diag(iatom)%mat(im1,im1,isppol,ispinor,ispinor)
518              end if
519 !            ! change it by: real(level_diag) with warning
520            end do
521          end do
522        end do
523 
524 !      ===================================
525 !      Compute possible occupations
526 !      ===================================
527 !      Value for nelec=0:
528        nconfig_nelec(0)=1
529        occ_level(0)%ocp(1,:)=0
530 !      Loop on possible occupation of levels with nelec
531        do nelec=1,nlevels ! number of electrons
532 !        write(message,'(2a,i3,a)') ch10," For number of electrons",  &
533 !        &       nelec," positions of electrons are:"
534 !        call wrtout(std_out,message,'COLL')
535 !        write(std_out,*) "nelec",nelec
536 !        write(std_out,*) "nlevels",nlevels
537          call combin(1,nconfig,nconfig_nelec,nelec,nlevels,occ_level,occup)
538          if(nconfig_nelec(nelec)/=nint(permutations(nlevels,nelec)/factorial(nelec))) then
539            message = " BUG in hubbard_one/combin"
540            MSG_BUG(message)
541          end if
542          occ_level(nelec)%ocp=zero
543          do iconfig=1,nconfig_nelec(nelec)
544            do ielec=1,nelec
545 !            occ_level%repart: gives the place of electron ielec for the configuration iconfig (among the config for the total number of electron nelec
546              occupied_level=occ_level(nelec)%repart(iconfig,ielec)
547 !            occ_level%ocp: gives if level occupied_level is occupied or not
548              occ_level(nelec)%ocp(iconfig,occupied_level)=1
549            end do
550          end do
551        end do
552 
553 !      ===================================
554 !      Print possible occupations
555 !      ===================================
556        if(prtopt>3) then
557          do nelec=0,nlevels ! number of electrons f
558            write(message,'(2a,i3,2a,i5,a)') ch10," For",nelec," electrons, ", &
559 &           "there are ",nconfig_nelec(nelec)," repartitions which are:"
560            call wrtout(std_out,message,'COLL')
561            do iconfig=1,nconfig_nelec(nelec)
562              write(message,'(40i4)') (occ_level(nelec)%ocp(iconfig,ilevel),ilevel=1,nlevels),&
563 &             (occ_level(nelec)%repart(iconfig,ielec),ielec=1,nelec)
564              call wrtout(std_out,message,'COLL')
565            end do
566          end do
567        end if
568 
569 !      ============================================
570 !      Compute energy for each of the occupations
571 !      ============================================
572        do nelec=0,nlevels !
573          e_nelec(nelec)%config=zero
574          do iconfig=1,nconfig_nelec(nelec)
575 !          First compute energy level contribution
576            do ielec=1,nelec
577              e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig) &
578 &             + elevels(occ_level(nelec)%repart(iconfig,ielec))
579            end do
580 !          write(std_out,*) "Nelec",nelec,"iconfig",iconfig,"eleve",e_nelec(nelec)%config(iconfig)
581 
582 !          Second: Compute interaction part
583 !          do ielec=1,nelec-1 ! compute interaction among the nelec electrons in the configuration iconfig
584 !          e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig)   &
585 !          &             + hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), &
586 !          &               occ_level(nelec)%repart(iconfig,ielec+1))
587 !          enddo
588            do ielec=1,nelec ! compute interaction among the nelec electrons in the configuration iconfig
589              do jelec=1,nelec
590                e_nelec(nelec)%config(iconfig)= e_nelec(nelec)%config(iconfig)   &
591 !              &               + hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), &
592 &               + udens_atoms(iatom)%value(occ_level(nelec)%repart(iconfig,ielec), &
593 &               occ_level(nelec)%repart(iconfig,jelec))/2.d0 ! udens(i,i)=0
594 !              write(std_out,*) ielec,occ_level(nelec)%repart(iconfig,ielec)
595 !              write(std_out,*) jelec,occ_level(nelec)%repart(iconfig,jelec)
596 !              write(std_out,*)hu(cryst_struc%typat(iatom))%udens(occ_level(nelec)%repart(iconfig,ielec), &
597 !              &                occ_level(nelec)%repart(iconfig,jelec))/2.d0
598              end do ! jelec
599            end do ! ielec
600 !          write(std_out,*) "Nelec",nelec,"iconfig",iconfig,"ecorr",e_nelec(nelec)%config(iconfig)
601 
602          end do ! iconfig
603          maxener(nelec)=maxval(-e_nelec(nelec)%config(:))
604          minener(nelec)=minval(-e_nelec(nelec)%config(:))
605        end do
606 !      write(std_out,*) "maxener", maxener(:)
607        emax=maxval(maxener(:))
608        emin=minval(minener(:))
609        eshift=zero
610        eshift=emax/two
611        eshift=emax-maxarg_exp/beta
612 !      eshift=emax
613 !      write(std_out,*)"emax",emax
614 !      write(std_out,*)"emin",emin
615 !      write(std_out,*)"eshift",eshift
616        write(message,'(a,3x,3a,3x,a)') ch10," Hubbard I: Energies as a", &
617 &       " function of number of electrons",ch10,&
618 &       "     Nelec     Min. Ene.       Max. Ener."
619        call wrtout(std_out,message,'COLL')
620        do nelec=0,nlevels
621          write(message,'(3x,a,i4,2f17.7)') "HI", nelec,&
622 &         minval(e_nelec(nelec)%config(:)),maxval(e_nelec(nelec)%config(:))
623          call wrtout(std_out,message,'COLL')
624        end do
625 
626 !      ===================================
627 !      Print possibles occupations
628 !      ===================================
629        if(prtopt>3) then
630          do nelec=0,nlevels ! number of electrons
631            write(message,'(2a,i3,2a,i5,3a)') ch10," For",nelec," electrons, ", &
632 &           "there are ",nconfig_nelec(nelec)," repartitions which are :", &
633 &           ch10,"Energy and Occupations"
634            call wrtout(std_out,message,'COLL')
635            do iconfig=1,nconfig_nelec(nelec)
636              write(message,'(f12.6,20i4)') e_nelec(nelec)%config(iconfig),&
637 &             (occ_level(nelec)%repart(iconfig,ielec),ielec=1,nelec)
638              call wrtout(std_out,message,'COLL')
639            end do
640          end do
641        end if
642 
643 !      sum_test=zero
644 !      do ielec=1,nelec+1
645 !      sum_test = sum_test + (occ_level(nelec)%repart(iconfig,ielec)  &
646 !      &              -occ_level(nelec)%repart(iconfig,ielec))
647 !      enddo
648 !      ===================================
649 !      Built transitions between configurations
650 !      ===================================
651        do nelec=0,nlevels-1
652          do iconfig=1,nconfig_nelec(nelec)
653            itrans=0 ! transition from iconfig
654            do jconfig=1, nconfig_nelec(nelec+1)
655              sum_test=0
656              do ilevel=1,nlevels
657 !              test if their is one electron added to the starting configuration
658                sum_test=sum_test + &
659 &               (occ_level(nelec+1)%ocp(jconfig,ilevel)- occ_level(nelec)%ocp(iconfig,ilevel))**2
660 !              save the level for the electron added
661                if(occ_level(nelec+1)%ocp(jconfig,ilevel)==1.and.occ_level(nelec)%ocp(iconfig,ilevel)==0) then
662                  m_temp=ilevel
663                end if
664              end do ! ilevel
665              if(sum_test==1) then
666                itrans=itrans+1
667                if(itrans>nlevels-nelec) then
668                  write(message,'(a,4i4)') "BUG: itrans is to big in hubbard_one",itrans,iconfig,jconfig,ilevel
669                  call wrtout(std_out,message,'COLL')
670                end if
671                occ_level(nelec)%transition(iconfig,itrans)=jconfig  ! jconfig=config(n+1) obtained after transition
672                occ_level(nelec)%transition_m(iconfig,itrans)=m_temp  !  level to fill to do the transition
673              end if
674            end do ! jconfig
675            if(prtopt>3) then
676              write(std_out,'(a,2i5,a,18i5)') "occ_level", nelec,&
677 &             iconfig,"  :",(occ_level(nelec)%transition(iconfig,itrans),itrans=1,nlevels-nelec)
678              write(std_out,'(a,2i5,a,18i5)') "electron added", nelec,iconfig,&
679 &             "  :",(occ_level(nelec)%transition_m(iconfig,itrans),itrans=1,nlevels-nelec)
680            end if
681          end do ! iconfig
682        end do ! nelec
683 
684 !      ===================================
685 !      Built Partition Function
686 !      ===================================
687        Z_part=czero
688 !      do nelec=1,nlevels-1
689        do nelec=0,nlevels
690          do iconfig=1,nconfig_nelec(nelec)
691            Ei_n    = e_nelec  (nelec  )%config(iconfig) + eshift
692            Z_part=Z_part+dexp(-Ei_n*beta)
693 !          write(std_out,*) "fonction de partition",nelec,iconfig, Z_part,Ei_n*beta,Ei_n,eshift
694          end do
695        end do
696 !      write(std_out,*) "Z_part",Z_part
697 
698 !      ===================================
699 !      Built Green Function
700 !      ===================================
701        ABI_ALLOCATE(green_temp,(green_hubbard%nw,nlevels))
702        ABI_ALLOCATE(green_temp_realw,(green_hubbard%nw,nlevels))
703 !      For each freq.
704 
705        green_temp=czero
706        green_temp_realw=czero
707        tmp=zero
708        do nelec=0,nlevels-1
709 !        write(std_out,*) "For nelec    =",nelec
710          do iconfig=1,nconfig_nelec(nelec)
711 !          write(std_out,*) "The config nb:",iconfig
712            do itrans=1,nlevels-nelec
713              jconfig = occ_level(nelec  )%transition(iconfig,itrans)
714              m_temp  = occ_level(nelec  )%transition_m(iconfig,itrans)
715              Ej_np1  = e_nelec  (nelec+1)%config(jconfig) + eshift
716              Ei_n    = e_nelec  (nelec  )%config(iconfig) + eshift
717 !            write(std_out,'(a,i4,a)') "Transition nb:",itrans,"involve"
718 !            write(std_out,'(a,i4,a)') "                        jconfig=",jconfig
719 !            write(std_out,'(a,i4,a)') "                        m_temp=",m_temp
720              do ifreq=1,green_hubbard%nw
721                if(green_hubbard%w_type=="imag") then
722                  omega_current=cmplx(zero,green_hubbard%omega(ifreq),kind=dp)
723                else if(green_hubbard%w_type=="real") then
724                  omega_current=cmplx(green_hubbard%omega(ifreq),zero,kind=dp)
725                end if
726                green_temp(ifreq,m_temp)=green_temp(ifreq,m_temp)+  &
727 &               (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta))/ &
728 &               ( omega_current +Ei_n-Ej_np1)
729                if(ifreq==1.and.m_temp==1) tmp=tmp+(dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta))
730 
731                green_temp_realw(ifreq,m_temp)=green_temp_realw(ifreq,m_temp)+  &
732 &               (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta))/ &
733 &               ( omega_current +Ei_n-Ej_np1)
734              end do
735 !            green_temp_realw(m_temp)=green_temp_realw(m_temp)+  &
736 !            &           (dexp(-Ej_np1*beta)+ dexp(-Ei_n*beta)) -> will give one at the end
737 !            write(std_out,*) "green",-Ej_np1*beta,-Ei_n*beta,dexp(-Ej_np1*beta),dexp(-Ei_n*beta)
738            end do
739          end do
740        end do
741 !      write(std_out,*) "tmp",tmp
742        ilevel=0
743        do ispinor=1,nspinor
744          do isppol=1,nsppol
745            do im=1,(2*lpawu+1)
746              ilevel=ilevel+1
747 !            write(std_out,'(16e15.6)') paw_dmft%omega_lo(ifreq),(real(green_temp_realw(ilevel)/Z_part),ilevel=1,nlevels)
748              do ifreq=1,green_hubbard%nw
749                green_hubbard%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=green_temp(ifreq,ilevel)/Z_part
750                green_hubbard_realw%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)=green_temp_realw(ifreq,ilevel)/Z_part
751              end do
752            end do
753          end do
754        end do
755 
756 !      End calculation for this frequency
757        ABI_DEALLOCATE(green_temp)
758        ABI_DEALLOCATE(green_temp_realw)
759 
760 !      ===================================
761 !      Deallocations
762 !      ===================================
763        do nelec=0,nlevels
764          ABI_DEALLOCATE(occ_level(nelec)%repart)
765          ABI_DEALLOCATE(occ_level(nelec)%ocp)
766          ABI_DEALLOCATE(occ_level(nelec)%transition)
767          ABI_DEALLOCATE(occ_level(nelec)%transition_m)
768          ABI_DEALLOCATE(e_nelec(nelec)%config)
769        end do
770        ABI_DATATYPE_DEALLOCATE(occ_level)
771        ABI_DEALLOCATE(occup)
772        ABI_DEALLOCATE(nconfig_nelec)
773        ABI_DATATYPE_DEALLOCATE(e_nelec)
774        ABI_DEALLOCATE(elevels)
775        ABI_DEALLOCATE(maxener)
776        ABI_DEALLOCATE(minener)
777      end if
778    end do
779    call destroy_green(green_hubbard_realw)
780 
781  end subroutine green_atomic_hubbard

m_hubbard_one/hubbard_one [ Functions ]

[ Top ] [ m_hubbard_one ] [ Functions ]

NAME

 hubbard_one

FUNCTION

 Solve the hubbard one approximation

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
  istep    =  step of iteration for LDA.
  lda_occup
  mpi_enreg=informations about MPI parallelization
  paw_dmft =  data for self-consistent LDA+DMFT calculations.
  pawang <type(pawang)>=paw angular mesh and related data

OUTPUT

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

NOTES

PARENTS

      impurity_solve,spectral_function

CHILDREN

      combin,destroy_green,init_green,wrtout

SOURCE

 81 subroutine hubbard_one(cryst_struc,green,hu,paw_dmft,pawang,pawprtvol,hdc,weiss)
 82 
 83  use defs_basis
 84  use m_errors
 85  use m_abicore
 86 
 87  use m_pawang, only : pawang_type
 88  use m_crystal, only : crystal_t
 89  use m_green, only : green_type,init_green,destroy_green
 90  use m_paw_dmft, only : paw_dmft_type
 91  use m_oper, only : oper_type,init_oper,destroy_oper,loc_oper,print_oper
 92  use m_matlu, only : matlu_type,sym_matlu, print_matlu, gather_matlu,&
 93 & diag_matlu,init_matlu,destroy_matlu,rotate_matlu,copy_matlu
 94  use m_hu, only : hu_type,rotatevee_hu
 95  use m_datafordmft, only : compute_levels
 96 
 97 !This section has been created automatically by the script Abilint (TD).
 98 !Do not modify the following lines by hand.
 99 #undef ABI_FUNC
100 #define ABI_FUNC 'hubbard_one'
101 !End of the abilint section
102 
103  implicit none
104 
105 !Arguments ------------------------------------
106 !scalars
107 ! type(pawang_type), intent(in) :: pawang
108  type(crystal_t),intent(in) :: cryst_struc
109  type(green_type), intent(inout) :: green
110  type(paw_dmft_type), intent(in)  :: paw_dmft
111  type(hu_type), intent(inout) :: hu(cryst_struc%ntypat)
112  type(pawang_type), intent(in) :: pawang
113  type(oper_type), intent(inout) :: hdc
114  integer, intent(in) :: pawprtvol
115  type(green_type), intent(inout) :: weiss
116 
117 !Local variables ------------------------------
118  type  :: level2_type
119   integer, pointer :: repart(:,:) => null()
120   integer, pointer :: ocp(:,:) => null()
121   integer, pointer :: transition(:,:) => null()
122   integer, pointer :: transition_m(:,:) => null()
123  end type level2_type
124  type  :: level1_type
125   real(dp), pointer :: config(:) => null()
126  end type level1_type
127 ! scalars
128  character(len=500) :: message
129  integer :: iatom,ifreq,im,im1,isppol,ispinor,ispinor1
130  integer :: lpawu,mbandc,natom,nkpt,nspinor,nsppol,nsppol_imp,tndim
131 ! complex(dpc) :: g,g0,w
132 ! arrays
133  complex(dp), allocatable :: Id(:,:,:,:)
134  type(coeff2c_type), allocatable :: eigvectmatlu(:,:)
135  type(coeff2_type), allocatable :: udens_atoms(:)
136  type(oper_type)  :: energy_level
137  type(green_type) :: green_hubbard
138  type(matlu_type), allocatable :: level_diag(:)
139  complex(dpc) :: omega_current
140 ! ************************************************************************
141  mbandc=paw_dmft%mbandc
142  nkpt=paw_dmft%nkpt
143  nsppol=paw_dmft%nsppol
144  natom=paw_dmft%natom
145  nspinor=paw_dmft%nspinor
146 
147 
148 !Initialise for compiler
149  omega_current=czero
150 
151 !======================================
152 !Allocations: levels and eigenvectors
153 !======================================
154  ABI_DATATYPE_ALLOCATE(level_diag,(natom))
155  ABI_DATATYPE_ALLOCATE(eigvectmatlu,(natom,nsppol))
156  ABI_DATATYPE_ALLOCATE(udens_atoms,(natom))
157  call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,level_diag)
158  do iatom=1,cryst_struc%natom
159    lpawu=paw_dmft%lpawu(iatom)
160    if(lpawu/=-1) then
161      tndim=nspinor*(2*lpawu+1)
162      do isppol=1,nsppol
163        ABI_ALLOCATE(eigvectmatlu(iatom,isppol)%value,(tndim,tndim))
164      end do
165      ABI_ALLOCATE(udens_atoms(iatom)%value,(2*(2*lpawu+1),2*(2*lpawu+1)))
166      level_diag(iatom)%mat=czero
167    end if
168  end do
169 
170  call init_oper(paw_dmft,energy_level,opt_ksloc=3)
171  call compute_levels(cryst_struc,energy_level,hdc,pawang,paw_dmft)
172 !!========================
173 !!Get KS eigenvalues
174 !!========================
175 ! call init_oper(paw_dmft,energy_level,opt_ksloc=3)
176 ! do iband=1,mbandc
177 !   do ikpt=1,nkpt
178 !     do isppol=1,nsppol
179 !!      Take \epsilon_{nks}
180 !!      ========================
181 !       energy_level%ks(isppol,ikpt,iband,iband)=paw_dmft%eigen_lda(isppol,ikpt,iband)
182 !     end do
183 !   end do
184 ! end do
185 !
186 !
187 !!======================================================================
188 !!Compute atomic levels from projection of \epsilon_{nks} and symetrize
189 !!======================================================================
190 ! call loc_oper(energy_level,paw_dmft,1)
191 ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels before sym and only LDA"
192 ! call wrtout(std_out,message,'COLL')
193 ! call print_matlu(energy_level%matlu,natom,1)
194 ! do iatom = 1 , natom
195 !   lpawu=paw_dmft%lpawu(iatom)
196 !   if(lpawu/=-1) then
197 !     do isppol=1,nsppol
198 !       do ispinor=1,nspinor
199 !         do im1=1,2*lpawu+1
200 !           energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)=&
201 !&           energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)&
202 !&           -hdc%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)-paw_dmft%fermie
203 !         end do
204 !       end do
205 !     end do
206 !!    write(std_out,*) "DC,fermie",hdc%matlu(iatom)%mat(1,1,1,1,1),paw_dmft%fermie
207 !   end if
208 ! end do ! natom
209 ! call sym_matlu(cryst_struc,energy_level%matlu,pawang)
210 !
211 ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels for Fermi Level=",paw_dmft%fermie
212 ! call wrtout(std_out,message,'COLL')
213 !!call print_oper(energy_level,1,paw_dmft,1)
214 ! call print_matlu(energy_level%matlu,natom,1)
215 
216 !========================
217 !Compute Weiss function
218 !========================
219  ABI_ALLOCATE(Id,(20,20,nspinor,nspinor))
220  do iatom = 1 , natom
221    lpawu=paw_dmft%lpawu(iatom)
222    if(lpawu/=-1) then
223      Id=czero
224      do im=1,2*lpawu+1
225        do ispinor=1,nspinor
226          Id(im,im,ispinor,ispinor)=cone
227        end do
228      end do ! ib
229      do ifreq=1,weiss%nw
230        if(weiss%w_type=="imag") then
231          omega_current=cmplx(zero,weiss%omega(ifreq),kind=dp)
232        else if(green%w_type=="real") then
233          omega_current=cmplx(weiss%omega(ifreq),zero,kind=dp)
234        end if
235        do im=1,2*lpawu+1
236          do im1=1,2*lpawu+1
237            do isppol=1,nsppol
238              do ispinor=1,nspinor
239                do ispinor1=1,nspinor
240                  weiss%oper(ifreq)%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=&
241 &                 ( omega_current*Id(im,im1,ispinor,ispinor1) - &
242 &                 energy_level%matlu(iatom)%mat(im,im1,isppol,ispinor,ispinor1))
243                end do ! ispinor1
244              end do ! ispinor
245            end do ! isppol
246          end do ! im1
247        end do ! im
248      end do ! ifreq
249    end if ! lpawu
250  end do ! natom
251  ABI_DEALLOCATE(Id)
252 
253 !=================================================================
254 !Diagonalizes atomic levels and keep eigenvectors in eigvectmatlu
255 !=================================================================
256 !if jpawu=0, rotatevee_hu will have no effect so it is not necessary to
257 !have a single rotation matrix for up and dn spins.
258 
259  if(hu(1)%jpawu_zero.and.nsppol==2) nsppol_imp=2
260  if(.not.hu(1)%jpawu_zero.or.nsppol/=2) nsppol_imp=1
261 !  Diagonalize energy levels
262  call diag_matlu(energy_level%matlu,level_diag,natom,&
263 & prtopt=pawprtvol,eigvectmatlu=eigvectmatlu,nsppol_imp=nsppol_imp,optreal=1)
264 
265 !  Use rotation matrix to rotate interaction
266  call rotatevee_hu(cryst_struc,hu,nspinor,nsppol,pawprtvol,eigvectmatlu,udens_atoms,1)
267 !write(std_out,*)"udens after rotatevee", udens_atoms(1)%value
268  write(message,'(a,2x,a,f13.5)') ch10,&
269 & " == Print Diagonalized Energy levels for Fermi Level=",paw_dmft%fermie
270  call wrtout(std_out,message,'COLL')
271  call print_matlu(level_diag,natom,1)
272 
273 !  Print out
274  if(nspinor==2) then
275    write(message,'(a,2x,a,f13.5)') ch10,&
276 &   " == Print weiss for small freq"
277    call wrtout(std_out,message,'COLL')
278    call print_matlu(weiss%oper(1)%matlu,natom,1)
279    write(message,'(a,2x,a,f13.5)') ch10,&
280 &   " == Print weiss for large freq"
281    call wrtout(std_out,message,'COLL')
282    call print_matlu(weiss%oper(weiss%nw)%matlu,natom,1)
283  end if
284 
285 !========================
286 !Compute Green function
287 !========================
288  call init_green(green_hubbard,paw_dmft,opt_oper_ksloc=2,wtype=green%w_type) ! initialize only matlu
289 !write(std_out,*)"udens", udens_atoms(1)%value
290 ! write(std_out,*)"levels",  level_diag(1)%mat
291  call green_atomic_hubbard(cryst_struc,green_hubbard,hu,level_diag,paw_dmft,udens_atoms)
292 !call rotate_matlu(energy_level%matlu,natom,pawprtvol=3)
293 !write(81,*) "I1",paw_dmft%omega_lo(1), real(green%oper(1)%matlu(1)%mat(1,1,1,1,1)),imag(green%oper(1)%matlu(1)%mat(1,1,1,1,1))
294 !========================================================================
295 !Rotate back Green function in the original basis before diagonalization
296 !========================================================================
297 !call print_matlu(level_diag,natom,1)
298 !test scall rotate_matlu(level_diag,eigvectmatlu,natom,3)
299 !todo_ab: add check here for back rotation
300 !call print_matlu(level_diag,natom,1)
301  write(message,'(2a,f13.5)') ch10," == Green function before rotation"
302  call wrtout(std_out,message,'COLL')
303  call print_matlu(green_hubbard%oper(1)%matlu,natom,1)
304  do ifreq=1,green_hubbard%nw
305    call rotate_matlu(green_hubbard%oper(ifreq)%matlu,eigvectmatlu,natom,3,0)
306    call copy_matlu(green_hubbard%oper(ifreq)%matlu,green%oper(ifreq)%matlu,natom)
307  end do
308  write(message,'(2a,f13.5)') ch10," == Green function after rotation"
309  call wrtout(std_out,message,'COLL')
310  call print_matlu(green%oper(1)%matlu,natom,1)
311  if(nspinor==2) then
312    write(message,'(a,2x,a,f13.5)') ch10,&
313 &   " == Print green for small freq"
314    call wrtout(std_out,message,'COLL')
315    call print_matlu(green%oper(1)%matlu,natom,1)
316    write(message,'(a,2x,a,f13.5)') ch10,&
317 &   " == Print green for large freq"
318    call wrtout(std_out,message,'COLL')
319    call print_matlu(green%oper(green%nw)%matlu,natom,1)
320  end if
321 !do ifreq=1,paw_dmft%dmft_nwlo
322 !g=green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
323 !g0=cone/weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)
324 !w=cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)
325 !write(160,*) paw_dmft%omega_lo(ifreq),real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)),imag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
326 !write(161,*) paw_dmft%omega_lo(ifreq),real(green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) ),imag(green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) )
327 !write(164,*) paw_dmft%omega_lo(ifreq),real(one/green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) ),imag(one/green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
328 !write(166,*) paw_dmft%omega_lo(ifreq),real(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)-one/green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) ),imag(weiss%oper(ifreq)%matlu(1)%mat(1,1,1,1,1)-one/green%oper(ifreq)%matlu(1)%mat(1,1,1,1,1))
329 !write(167,*) paw_dmft%omega_lo(ifreq),real((g-g0)/(g0*g)),imag((g-g0)/(g0*g))
330 !write(168,*) paw_dmft%omega_lo(ifreq),real(1/g-w),imag(1/g-w)
331 !write(169,*) paw_dmft%omega_lo(ifreq),real(1/g0-w),imag(1/g0-w)
332 !write(170,*) paw_dmft%omega_lo(ifreq),w
333 !write(171,*) paw_dmft%omega_lo(ifreq),real(1/g),imag(1/g)
334 !write(172,*) paw_dmft%omega_lo(ifreq),real(w),imag(w)
335 
336 !! voir si en faisant GG0/(G-G0) cela reduit l'erreur
337 !enddo
338 !call abi_abort('COLL')
339 
340 
341 !write(message,'(2a,f13.5)') ch10," == Print Energy levels after diagonalisation"
342 !call wrtout(std_out,message,'COLL')
343 !call print_matlu(energy_level%matlu,natom,1)
344 
345 !======================================
346 !Deallocations and destroys
347 !======================================
348  call destroy_green(green_hubbard)
349  call destroy_oper(energy_level)
350  call destroy_matlu(level_diag,natom)
351  ABI_DATATYPE_DEALLOCATE(level_diag)
352  do iatom=1,cryst_struc%natom
353    lpawu=paw_dmft%lpawu(iatom)
354    if(lpawu/=-1) then
355      ABI_DEALLOCATE(udens_atoms(iatom)%value)
356      do isppol=1,nsppol
357        ABI_DEALLOCATE(eigvectmatlu(iatom,isppol)%value)
358      end do
359    end if
360  end do
361  ABI_DATATYPE_DEALLOCATE(eigvectmatlu)
362  ABI_DATATYPE_DEALLOCATE(udens_atoms)