diff --git a/src/ED_AUX_FUNX.f90 b/src/ED_AUX_FUNX.f90 index cf94e6a3..3315d7df 100644 --- a/src/ED_AUX_FUNX.f90 +++ b/src/ED_AUX_FUNX.f90 @@ -12,8 +12,11 @@ MODULE ED_AUX_FUNX !> ED SET HLOC interface ed_set_Hloc - module procedure :: ed_set_Hloc_single - module procedure :: ed_set_Hloc_lattice + module procedure :: ed_set_Hloc_single_N2 + module procedure :: ed_set_Hloc_single_N4 + module procedure :: ed_set_Hloc_lattice_N2 + module procedure :: ed_set_Hloc_lattice_N3 + module procedure :: ed_set_Hloc_lattice_N5 end interface ed_set_Hloc interface lso2nnn_reshape @@ -164,9 +167,9 @@ end subroutine ed_set_suffix_c !+------------------------------------------------------------------+ !PURPOSE : Setup Himpurity, the local part of the non-interacting Hamiltonian - !+------------------------------------------------------------------+ - subroutine ed_set_Hloc_single(Hloc) - complex(8),dimension(..),intent(in) :: Hloc + !+------------------------------------------------------------------+ + subroutine ed_set_Hloc_single_N2(Hloc) + complex(8),dimension(:,:),intent(in) :: Hloc #ifdef _DEBUG write(Logfile,"(A)")"DEBUG ed_set_Hloc: set impHloc" #endif @@ -174,22 +177,28 @@ subroutine ed_set_Hloc_single(Hloc) if(allocated(impHloc))deallocate(impHloc) allocate(impHloc(Nspin,Nspin,Norb,Norb));impHloc=zero ! - select rank(Hloc) - rank default;stop "ED_SET_HLOC ERROR: Hloc has a wrong rank. Accepted: [Nso,Nso] or [Nspin,Nspin,Norb,Norb]" - rank (2) !Hloc[Nso,Nso] call assert_shape(Hloc,[Nspin*Norb,Nspin*Norb],"ed_set_Hloc","Hloc") impHloc = so2nn_reshape(Hloc(1:Nspin*Norb,1:Nspin*Norb),Nspin,Norb) - rank (4) !Hloc[Nspin,Nspin,Norb,Norb] + if(ed_verbose>2)call print_hloc(impHloc) + end subroutine ed_set_Hloc_single_N2 + + subroutine ed_set_Hloc_single_N4(Hloc) + complex(8),dimension(:,:,:,:),intent(in) :: Hloc +#ifdef _DEBUG + write(Logfile,"(A)")"DEBUG ed_set_Hloc: set impHloc" +#endif + ! + if(allocated(impHloc))deallocate(impHloc) + allocate(impHloc(Nspin,Nspin,Norb,Norb));impHloc=zero + ! call assert_shape(Hloc,[Nspin,Nspin,Norb,Norb],"ed_set_Hloc","Hloc") impHloc = Hloc - end select if(ed_verbose>2)call print_hloc(impHloc) - end subroutine ed_set_Hloc_single - + end subroutine ed_set_Hloc_single_N4 - subroutine ed_set_Hloc_lattice(Hloc,Nlat) - complex(8),dimension(..),intent(in) :: Hloc - integer :: Nlat,ilat + subroutine ed_set_Hloc_lattice_N2(Hloc,Nlat) + complex(8),dimension(:,:),intent(in) :: Hloc + integer :: Nlat,ilat #ifdef _DEBUG write(Logfile,"(A)")"DEBUG ed_set_Hloc: set impHloc" #endif @@ -197,27 +206,41 @@ subroutine ed_set_Hloc_lattice(Hloc,Nlat) if(allocated(Hloc_ineq))deallocate(Hloc_ineq) allocate(Hloc_ineq(Nlat,Nspin,Nspin,Norb,Norb));Hloc_ineq=zero ! - select rank(Hloc) - rank default; - stop "ED_SET_HLOC ERROR: Hloc has a wrong rank. [Nlso,Nlso],[Nlat,Nso,Nso],[Nlat,Nspin,Nspin,Norb,Norb]" - ! - rank (2) call assert_shape(Hloc,[Nlat*Nspin*Norb,Nlat*Nspin*Norb],'ed_set_Hloc','Hloc') Hloc_ineq = lso2nnn_reshape(Hloc(1:Nlat*Nspin*Norb,1:Nlat*Nspin*Norb),Nlat,Nspin,Norb) + end subroutine ed_set_Hloc_lattice_N2 + + + subroutine ed_set_Hloc_lattice_N3(Hloc,Nlat) + complex(8),dimension(:,:,:),intent(in) :: Hloc + integer :: Nlat,ilat +#ifdef _DEBUG + write(Logfile,"(A)")"DEBUG ed_set_Hloc: set impHloc" +#endif + ! + if(allocated(Hloc_ineq))deallocate(Hloc_ineq) + allocate(Hloc_ineq(Nlat,Nspin,Nspin,Norb,Norb));Hloc_ineq=zero ! - rank (3) call assert_shape(Hloc,[Nlat,Nspin*Norb,Nspin*Norb],'ed_set_Hloc','Hloc') do ilat=1,Nlat Hloc_ineq(ilat,:,:,:,:) = so2nn_reshape(Hloc(ilat,1:Nspin*Norb,1:Nspin*Norb),Nspin,Norb) enddo ! - rank (5) + end subroutine ed_set_Hloc_lattice_N3 + + subroutine ed_set_Hloc_lattice_N5(Hloc,Nlat) + complex(8),dimension(:,:,:,:,:),intent(in) :: Hloc + integer :: Nlat,ilat +#ifdef _DEBUG + write(Logfile,"(A)")"DEBUG ed_set_Hloc: set impHloc" +#endif + ! + if(allocated(Hloc_ineq))deallocate(Hloc_ineq) + allocate(Hloc_ineq(Nlat,Nspin,Nspin,Norb,Norb));Hloc_ineq=zero + ! call assert_shape(Hloc,[Nlat,Nspin,Nspin,Norb,Norb],'ed_set_Hloc','Hloc') Hloc_ineq = Hloc - end select - end subroutine ed_set_Hloc_lattice - - + end subroutine ed_set_Hloc_lattice_N5 diff --git a/src/ED_BATH/ED_FIT_BATH.f90 b/src/ED_BATH/ED_FIT_BATH.f90 index 1f59fe4b..c2ef993b 100644 --- a/src/ED_BATH/ED_FIT_BATH.f90 +++ b/src/ED_BATH/ED_FIT_BATH.f90 @@ -25,10 +25,16 @@ MODULE ED_BATH_FIT private interface ed_chi2_fitgf - module procedure chi2_fitgf_single_normal - module procedure chi2_fitgf_single_superc - module procedure chi2_fitgf_lattice_normal - module procedure chi2_fitgf_lattice_superc + module procedure chi2_fitgf_single_normal_n3 + module procedure chi2_fitgf_single_normal_n5 + module procedure chi2_fitgf_single_superc_n3 + module procedure chi2_fitgf_single_superc_n5 + module procedure chi2_fitgf_lattice_normal_n3 + module procedure chi2_fitgf_lattice_normal_n4 + module procedure chi2_fitgf_lattice_normal_n6 + module procedure chi2_fitgf_lattice_superc_n3 + module procedure chi2_fitgf_lattice_superc_n4 + module procedure chi2_fitgf_lattice_superc_n6 end interface ed_chi2_fitgf public :: ed_chi2_fitgf @@ -43,15 +49,15 @@ MODULE ED_BATH_FIT ! - CHI2_FITGF_GENERIC_NORMAL interface for the normal case ! * CHI2_FITGF_GENERIC_NORMAL_NOSPIN interface to fixed spin input !+----------------------------------------------------------------------+ - subroutine chi2_fitgf_single_normal(g,bath,ispin,iorb,fmpi) - complex(8),dimension(..) :: g - real(8),dimension(:) :: bath - integer,optional :: ispin,iorb - logical,optional :: fmpi + subroutine chi2_fitgf_single_normal_n3(g,bath,ispin,iorb,fmpi) + complex(8),dimension(:,:,:) :: g + real(8),dimension(:) :: bath + integer,optional :: ispin,iorb + logical,optional :: fmpi ! complex(8),dimension(Nspin,Nspin,Norb,Norb,Lmats) :: fg - integer :: ispin_,Liw - logical :: fmpi_ + integer :: ispin_,Liw + logical :: fmpi_ ! write(Logfile,"(A)")"" #ifdef _DEBUG @@ -76,16 +82,107 @@ subroutine chi2_fitgf_single_normal(g,bath,ispin,iorb,fmpi) ! if(MpiMaster)then ! - select rank(g) - rank default - stop "chi2_fitgf_generic_normal error: input G has a wrong rank: accepted 3 or 5" - rank (3) call assert_shape(g,[Nspin*Norb,Nspin*Norb,Lmats],"chi2_fitgf_generic_normal","g") fg = so2nn_reshape(g(1:Nspin*Norb,1:Nspin*Norb,1:Lmats),Nspin,Norb,Lmats) - rank (5) + ! + select case(bath_type) + case default + select case(ed_mode) + case ("normal") + if(present(iorb))then + call chi2_fitgf_normal_normal(fg(ispin_,ispin_,:,:,:),bath,ispin_,iorb) + else + call chi2_fitgf_normal_normal(fg(ispin_,ispin_,:,:,:),bath,ispin_) + endif + case ("nonsu2") + call chi2_fitgf_normal_nonsu2(fg(:,:,:,:,:),bath) + case default + stop "chi2_fitgf ERROR: ed_mode!=normal/nonsu2 but only NORMAL component is provided" + end select + ! + case ("hybrid") + select case(ed_mode) + case ("normal") + call chi2_fitgf_hybrid_normal(fg(ispin_,ispin_,:,:,:),bath,ispin_) + case ("nonsu2") + call chi2_fitgf_hybrid_nonsu2(fg(:,:,:,:,:),bath) + case default + stop "chi2_fitgf ERROR: ed_mode!=normal/nonsu2 but only NORMAL component is provided" + end select + ! + case ("replica") + select case(ed_mode) + case ("normal","nonsu2") + call chi2_fitgf_replica(fg,bath) + case default + stop "chi2_fitgf ERROR: ed_mode!=normal/nonsu2 but only NORMAL component is provided" + end select + case ("general") + select case(ed_mode) + case ("normal","nonsu2") + call chi2_fitgf_general(fg,bath) + case default + stop "chi2_fitgf ERROR: ed_mode!=normal/nonsu2 but only NORMAL component is provided" + end select + end select + endif + ! +#ifdef _MPI + if(MpiStatus)then + call Bcast_MPI(MpiComm,bath) + if(.not.MpiMaster)write(LOGfile,"(A)")"Bath received from master node" + endif +#endif + ! + !set trim_state_list to true after the first fit has been done: this + !marks the ends of the cycle of the 1st DMFT loop. + trim_state_list=.true. + !DELETE THE LOCAL MPI COMMUNICATOR: +#ifdef _MPI + if(check_MPI().AND.fmpi_)call ed_del_MpiComm() +#endif +#ifdef _DEBUG + write(Logfile,"(A)")"" +#endif + end subroutine chi2_fitgf_single_normal_n3 + + + + subroutine chi2_fitgf_single_normal_n5(g,bath,ispin,iorb,fmpi) + complex(8),dimension(:,:,:,:,:) :: g + real(8),dimension(:) :: bath + integer,optional :: ispin,iorb + logical,optional :: fmpi + ! + complex(8),dimension(Nspin,Nspin,Norb,Norb,Lmats) :: fg + integer :: ispin_,Liw + logical :: fmpi_ + ! + write(Logfile,"(A)")"" +#ifdef _DEBUG + write(Logfile,"(A)")"DEBUG chi2_fitgf_generic_normal_mpi: Start Chi**2 fit" +#endif + ! + ispin_=1;if(present(ispin))ispin_=ispin + fmpi_=.true.;if(present(fmpi))fmpi_=fmpi + ! +#ifdef _MPI + if(check_MPI().AND.fmpi_)call ed_set_MpiComm() +#endif + ! + select case(cg_method) + case default + stop "ED Error: cg_method > 1" + case (0) + if(ed_verbose>2)write(LOGfile,"(A,I1,A,A)")"Chi^2 fit with CG-nr and CG-weight: ",cg_weight," on: ",cg_scheme + case (1) + if(ed_verbose>2)write(LOGfile,"(A,I1,A,A)")"Chi^2 fit with CG-minimize and CG-weight: ",cg_weight," on: ",cg_scheme + end select + ! + if(MpiMaster)then + ! call assert_shape(g,[Nspin,Nspin,Norb,Norb,Lmats],"chi2_fitgf_generic_normal","g") fg = g(1:Nspin,1:Nspin,1:Norb,1:Norb,1:Lmats) - end select ! select case(bath_type) case default @@ -146,12 +243,16 @@ subroutine chi2_fitgf_single_normal(g,bath,ispin,iorb,fmpi) #ifdef _DEBUG write(Logfile,"(A)")"" #endif - end subroutine chi2_fitgf_single_normal + end subroutine chi2_fitgf_single_normal_n5 + + + + + + - - subroutine chi2_fitgf_single_superc(g,f,bath,ispin,iorb,fmpi) - complex(8),dimension(..) :: g - complex(8),dimension(..) :: f + subroutine chi2_fitgf_single_superc_n3(g,f,bath,ispin,iorb,fmpi) + complex(8),dimension(:,:,:) :: g,f real(8),dimension(:) :: bath integer,optional :: ispin,iorb logical,optional :: fmpi @@ -184,27 +285,108 @@ subroutine chi2_fitgf_single_superc(g,f,bath,ispin,iorb,fmpi) ! if(MpiMaster)then ! - select rank(g) - rank default - stop "chi2_fitgf_generic_superc error: input G has a wrong rank: accepted 3 or 5" - rank (3) call assert_shape(g,[Nspin*Norb,Nspin*Norb,Lmats],"chi2_fitgf_generic_superc","g") fg(1,:,:,:,:,:) = so2nn_reshape(g(1:Nspin*Norb,1:Nspin*Norb,1:Lmats),Nspin,Norb,Lmats) - rank (5) - call assert_shape(g,[Nspin,Nspin,Norb,Norb,Lmats],"chi2_fitgf_generic_superc","g") - fg(1,:,:,:,:,:) = g(1:Nspin,1:Nspin,1:Norb,1:Norb,1:Lmats) - end select - ! - select rank(f) - rank default - stop "chi2_fitgf_generic_superc error: input F has a wrong rank: accepted 3 or 5" - rank (3) call assert_shape(f,[Nspin*Norb,Nspin*Norb,Lmats],"chi2_fitgf_generic_superc","f") fg(2,:,:,:,:,:) = so2nn_reshape(f(1:Nspin*Norb,1:Nspin*Norb,1:Lmats),Nspin,Norb,Lmats) - rank (5) + ! + select case(bath_type) + case default + select case(ed_mode) + case ("superc") + if(present(iorb))then + call chi2_fitgf_normal_superc(fg(:,ispin_,ispin_,:,:,:),bath,ispin_,iorb) + else + call chi2_fitgf_normal_superc(fg(:,ispin_,ispin_,:,:,:),bath,ispin_) + endif + case default + write(LOGfile,"(A)") "chi2_fitgf WARNING: ed_mode=normal/nonsu2 but NORMAL & ANOMAL components provided." + call chi2_fitgf_normal_normal(fg(1,ispin_,ispin_,:,:,:),bath,ispin_) + end select + case ("hybrid") + select case(ed_mode) + case ("superc") + call chi2_fitgf_hybrid_superc(fg(:,ispin_,ispin_,:,:,:),bath,ispin_) + case default + write(LOGfile,"(A)") "chi2_fitgf WARNING: ed_mode=normal/nonsu2 but NORMAL & ANOMAL components provided." + call chi2_fitgf_hybrid_normal(fg(1,ispin_,ispin_,:,:,:),bath,ispin_) + end select + case ("replica") + select case(ed_mode) + case ("superc") + call chi2_fitgf_replica_superc(fg(:,:,:,:,:,:),bath) + case default + write(LOGfile,"(A)") "chi2_fitgf WARNING: ed_mode=normal/nonsu2 but NORMAL & ANOMAL components provided." + call chi2_fitgf_replica(fg(1,:,:,:,:,:),bath) + end select + case ("general") + select case(ed_mode) + case ("superc") + call chi2_fitgf_general_superc(fg(:,:,:,:,:,:),bath) + case default + write(LOGfile,"(A)") "chi2_fitgf WARNING: ed_mode=normal/nonsu2 but NORMAL & ANOMAL components provided." + call chi2_fitgf_general(fg(1,:,:,:,:,:),bath) + end select + end select + endif + ! +#ifdef _MPI + if(MpiStatus)then + call Bcast_MPI(MpiComm,bath) + if(.not.MpiMaster)write(LOGfile,"(A)")"Bath received from master node" + endif +#endif + ! + !set trim_state_list to true after the first fit has been done: this + !marks the ends of the cycle of the 1st DMFT loop. + trim_state_list=.true. +#ifdef _MPI + if(check_MPI().AND.fmpi_)call ed_del_MpiComm() +#endif +#ifdef _DEBUG + write(Logfile,"(A)")"" +#endif + end subroutine chi2_fitgf_single_superc_n3 + + + subroutine chi2_fitgf_single_superc_n5(g,f,bath,ispin,iorb,fmpi) + complex(8),dimension(:,:,:,:,:) :: g,f + real(8),dimension(:) :: bath + integer,optional :: ispin,iorb + logical,optional :: fmpi + ! + complex(8),dimension(2,Nspin,Nspin,Norb,Norb,Lmats) :: fg + integer :: ispin_ + logical :: fmpi_ + ! + write(Logfile,"(A)")"" +#ifdef _DEBUG + write(Logfile,"(A)")"DEBUG chi2_fitgf_generic_superc: Start Chi**2 fit" +#endif + ! + ! + ispin_=1;if(present(ispin))ispin_=ispin + fmpi_=.true.;if(present(fmpi))fmpi_=fmpi + ! +#ifdef _MPI + if(check_MPI().AND.fmpi_)call ed_set_MpiComm() +#endif + ! + select case(cg_method) + case default + stop "ED Error: cg_method > 1" + case (0) + if(ed_verbose>2)write(LOGfile,"(A,I1,A,A)")"master: Chi^2 fit with CG-nr and CG-weight: ",cg_weight," on: ",cg_scheme + case (1) + if(ed_verbose>2)write(LOGfile,"(A,I1,A,A)")"master: Chi^2 fit with CG-minimize and CG-weight: ",cg_weight," on: ",cg_scheme + end select + ! + if(MpiMaster)then + ! + call assert_shape(g,[Nspin,Nspin,Norb,Norb,Lmats],"chi2_fitgf_generic_superc","g") + fg(1,:,:,:,:,:) = g(1:Nspin,1:Nspin,1:Norb,1:Norb,1:Lmats) call assert_shape(f,[Nspin,Nspin,Norb,Norb,Lmats],"chi2_fitgf_generic_superc","f") fg(2,:,:,:,:,:) = f(1:Nspin,1:Nspin,1:Norb,1:Norb,1:Lmats) - end select ! select case(bath_type) case default @@ -262,33 +444,36 @@ subroutine chi2_fitgf_single_superc(g,f,bath,ispin,iorb,fmpi) #ifdef _DEBUG write(Logfile,"(A)")"" #endif - end subroutine chi2_fitgf_single_superc - + end subroutine chi2_fitgf_single_superc_n5 + !################################################################## + !################################################################## + !################################################################## + !################################################################## - subroutine chi2_fitgf_lattice_normal(g,bath,ispin) - real(8),intent(inout) :: bath(:,:) - complex(8),dimension(..) :: g - integer,optional :: ispin + subroutine chi2_fitgf_lattice_normal_n3(g,bath,ispin) + real(8),intent(inout) :: bath(:,:) + complex(8),dimension(:,:,:) :: g + integer,optional :: ispin ! - complex(8) :: fg(size(bath,1),Nspin,Nspin,Norb,Norb,Lmats) - real(8) :: bath_tmp(size(bath,1),size(bath,2)) - integer :: i,ispin_ - integer :: ilat - integer :: iorb,is,io - integer :: jorb,js,jo - integer :: Nsites - logical :: check_dim - character(len=5) :: tmp_suffix - integer :: MPI_ID=0 - integer :: MPI_SIZE=1 - logical :: MPI_MASTER=.true. + complex(8) :: fg(size(bath,1),Nspin,Nspin,Norb,Norb,Lmats) + real(8) :: bath_tmp(size(bath,1),size(bath,2)) + integer :: i,ispin_ + integer :: ilat + integer :: iorb,is,io + integer :: jorb,js,jo + integer :: Nsites + logical :: check_dim + character(len=5) :: tmp_suffix + integer :: MPI_ID=0 + integer :: MPI_SIZE=1 + logical :: MPI_MASTER=.true. ! write(Logfile,"(A)")"" ! @@ -313,27 +498,159 @@ subroutine chi2_fitgf_lattice_normal(g,bath,ispin) ! fg = zero ! - select rank(g) - rank default - stop "chi2_fitgf_generic_normal error: input G has a wrong rank: accepted 3 or 5" - ! - rank (3) if(ed_verbose>1)write(Logfile,"(A)")"Chi**2 get G with rank:"//str(rank(g)) call assert_shape(g,[Nsites*Nspin*Norb,Nsites*Nspin*Norb,Lmats],'chi2_fitgf_generic_normal','g') fg = lso2nnn_reshape(g(1:Nsites*Nspin*Norb,1:Nsites*Nspin*Norb,1:Lmats),Nsites,Nspin,Norb,Lmats) ! - rank (4) + ! + bath_tmp=0d0 + do ilat = 1+MPI_ID,Nsites,MPI_SIZE + if(ed_verbose>1)write(Logfile,"(A)")"Start Chi**2 fit for site:"//str(ilat) + bath_tmp(ilat,:)=bath(ilat,:) + call ed_set_Hloc(Hloc_ineq(ilat,:,:,:,:)) + ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) + ! + if(present(ispin))then + call ed_chi2_fitgf(fg(ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin=ispin_,fmpi=.false.) + else + call ed_chi2_fitgf(fg(ilat,:,:,:,:,:),bath_tmp(ilat,:),fmpi=.false.) + end if + end do + ! +#ifdef _MPI + if(check_MPI())then + bath=0d0 + call AllReduce_MPI(MPI_COMM_WORLD,bath_tmp,bath) + else + bath = bath_tmp + endif +#else + bath = bath_tmp +#endif + ! + call ed_reset_suffix + end subroutine chi2_fitgf_lattice_normal_n3 + + subroutine chi2_fitgf_lattice_normal_n4(g,bath,ispin) + real(8),intent(inout) :: bath(:,:) + complex(8),dimension(:,:,:,:) :: g + integer,optional :: ispin + ! + complex(8) :: fg(size(bath,1),Nspin,Nspin,Norb,Norb,Lmats) + real(8) :: bath_tmp(size(bath,1),size(bath,2)) + integer :: i,ispin_ + integer :: ilat + integer :: iorb,is,io + integer :: jorb,js,jo + integer :: Nsites + logical :: check_dim + character(len=5) :: tmp_suffix + integer :: MPI_ID=0 + integer :: MPI_SIZE=1 + logical :: MPI_MASTER=.true. + ! + write(Logfile,"(A)")"" + ! + ispin_=1;if(present(ispin))ispin_=ispin + if(ispin_>Nspin)stop "ed_fit_bath_sites error: required spin index > Nspin" + ! +#ifdef _MPI + if(check_MPI())then + MPI_ID = get_Rank_MPI() + MPI_SIZE = get_Size_MPI() + MPI_MASTER = get_Master_MPI() + endif +#endif + ! + ! Check dimensions ! + Nsites=size(bath,1) + ! + do ilat=1+MPI_ID,Nsites,MPI_SIZE + check_dim = check_bath_dimension(bath(ilat,:)) + if(.not.check_dim) stop "init_lattice_bath: wrong bath size dimension 1 or 2 " + end do + ! + fg = zero + ! if(ed_verbose>1)write(Logfile,"(A)")"Chi**2 get G with rank:"//str(rank(g)) call assert_shape(g,[Nsites,Nspin*Norb,Nspin*Norb,Lmats],'chi2_fitgf_generic_normal','g') do ilat=1,Nsites fg(ilat,:,:,:,:,:) = so2nn_reshape(g(ilat,1:Nspin*Norb,1:Nspin*Norb,1:Lmats),Nspin,Norb,Lmats) enddo ! - rank (6) + ! + bath_tmp=0d0 + do ilat = 1+MPI_ID,Nsites,MPI_SIZE + if(ed_verbose>1)write(Logfile,"(A)")"Start Chi**2 fit for site:"//str(ilat) + bath_tmp(ilat,:)=bath(ilat,:) + call ed_set_Hloc(Hloc_ineq(ilat,:,:,:,:)) + ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) + ! + if(present(ispin))then + call ed_chi2_fitgf(fg(ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin=ispin_,fmpi=.false.) + else + call ed_chi2_fitgf(fg(ilat,:,:,:,:,:),bath_tmp(ilat,:),fmpi=.false.) + end if + end do + ! +#ifdef _MPI + if(check_MPI())then + bath=0d0 + call AllReduce_MPI(MPI_COMM_WORLD,bath_tmp,bath) + else + bath = bath_tmp + endif +#else + bath = bath_tmp +#endif + ! + call ed_reset_suffix + end subroutine chi2_fitgf_lattice_normal_n4 + + subroutine chi2_fitgf_lattice_normal_n6(g,bath,ispin) + real(8),intent(inout) :: bath(:,:) + complex(8),dimension(:,:,:,:,:,:) :: g + integer,optional :: ispin + ! + complex(8) :: fg(size(bath,1),Nspin,Nspin,Norb,Norb,Lmats) + real(8) :: bath_tmp(size(bath,1),size(bath,2)) + integer :: i,ispin_ + integer :: ilat + integer :: iorb,is,io + integer :: jorb,js,jo + integer :: Nsites + logical :: check_dim + character(len=5) :: tmp_suffix + integer :: MPI_ID=0 + integer :: MPI_SIZE=1 + logical :: MPI_MASTER=.true. + ! + write(Logfile,"(A)")"" + ! + ispin_=1;if(present(ispin))ispin_=ispin + if(ispin_>Nspin)stop "ed_fit_bath_sites error: required spin index > Nspin" + ! +#ifdef _MPI + if(check_MPI())then + MPI_ID = get_Rank_MPI() + MPI_SIZE = get_Size_MPI() + MPI_MASTER = get_Master_MPI() + endif +#endif + ! + ! Check dimensions ! + Nsites=size(bath,1) + ! + do ilat=1+MPI_ID,Nsites,MPI_SIZE + check_dim = check_bath_dimension(bath(ilat,:)) + if(.not.check_dim) stop "init_lattice_bath: wrong bath size dimension 1 or 2 " + end do + ! + fg = zero + ! if(ed_verbose>1)write(Logfile,"(A)")"Chi**2 get G with rank:"//str(rank(g)) call assert_shape(g,[Nsites,Nspin,Nspin,Norb,Norb,Lmats],'chi2_fitgf_generic_normal','g') fg = g(1:Nsites,1:Nspin,1:Nspin,1:Norb,1:Norb,1:Lmats) - end select ! bath_tmp=0d0 do ilat = 1+MPI_ID,Nsites,MPI_SIZE @@ -361,26 +678,26 @@ subroutine chi2_fitgf_lattice_normal(g,bath,ispin) #endif ! call ed_reset_suffix - end subroutine chi2_fitgf_lattice_normal - - - - - - subroutine chi2_fitgf_lattice_superc(g,f,bath,ispin) - real(8),intent(inout) :: bath(:,:) - complex(8),dimension(..) :: g,f - integer,optional :: ispin - ! - complex(8) :: fg(2,size(bath,1),Nspin,Nspin,Norb,Norb,Lmats) - real(8) :: bath_tmp(size(bath,1),size(bath,2)) - integer :: ilat,i,iorb,ispin_ - integer :: Nsites - logical :: check_dim - character(len=5) :: tmp_suffix - integer :: MPI_ID=0 - integer :: MPI_SIZE=1 - logical :: MPI_MASTER=.true. + end subroutine chi2_fitgf_lattice_normal_n6 + + + + + + subroutine chi2_fitgf_lattice_superc_n3(g,f,bath,ispin) + real(8),intent(inout) :: bath(:,:) + complex(8),dimension(:,:,:) :: g,f + integer,optional :: ispin + ! + complex(8) :: fg(2,size(bath,1),Nspin,Nspin,Norb,Norb,Lmats) + real(8) :: bath_tmp(size(bath,1),size(bath,2)) + integer :: ilat,i,iorb,ispin_ + integer :: Nsites + logical :: check_dim + character(len=5) :: tmp_suffix + integer :: MPI_ID=0 + integer :: MPI_SIZE=1 + logical :: MPI_MASTER=.true. ! write(Logfile,"(A)")"" ! @@ -403,41 +720,172 @@ subroutine chi2_fitgf_lattice_superc(g,f,bath,ispin) if(.not.check_dim) stop "init_lattice_bath: wrong bath size dimension 1 or 2 " end do ! - select rank(g) - rank default - stop "chi2_fitgf_generic_superc error: input G has a wrong rank: accepted 3 or 5" - rank (3) call assert_shape(g,[Nsites*Nspin*Norb,Nsites*Nspin*Norb,Lmats],'chi2_fitgf_generic_superc','g') fg(1,:,:,:,:,:,:) = lso2nnn_reshape(g(1:Nsites*Nspin*Norb,1:Nsites*Nspin*Norb,1:Lmats),Nsites,Nspin,Norb,Lmats) ! - rank (4) + call assert_shape(f,[Nsites*Nspin*Norb,Nsites*Nspin*Norb,Lmats],'chi2_fitgf_generic_superc','f') + fg(2,:,:,:,:,:,:) = lso2nnn_reshape(f(1:Nsites*Nspin*Norb,1:Nsites*Nspin*Norb,1:Lmats),Nsites,Nspin,Norb,Lmats) + ! + ! + bath_tmp=0.d0 + ! + do ilat= 1 + MPI_ID, Nsites, MPI_SIZE + write(Logfile,"(A)")"ed_fit_bath_sites_superc: Start Chi**2 fit for site:"//str(ilat) + bath_tmp(ilat,:) = bath(ilat,:) + call ed_set_Hloc(Hloc_ineq(ilat,:,:,:,:)) + ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) + ! + if(present(ispin))then + call ed_chi2_fitgf(fg(1,ilat,:,:,:,:,:),fg(2,ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin_,fmpi=.false.) + else + do ispin_=1,Nspin + call ed_chi2_fitgf(fg(1,ilat,:,:,:,:,:),fg(2,ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin_,fmpi=.false.) + enddo + endif + end do + ! +#ifdef _MPI + if(check_MPI())then + bath=0d0 + call AllReduce_MPI(MPI_COMM_WORLD,bath_tmp,bath) + else + bath = bath_tmp + endif +#else + bath = bath_tmp +#endif + ! + call ed_reset_suffix + end subroutine chi2_fitgf_lattice_superc_n3 + + + + + + + subroutine chi2_fitgf_lattice_superc_n4(g,f,bath,ispin) + real(8),intent(inout) :: bath(:,:) + complex(8),dimension(:,:,:,:) :: g,f + integer,optional :: ispin + ! + complex(8) :: fg(2,size(bath,1),Nspin,Nspin,Norb,Norb,Lmats) + real(8) :: bath_tmp(size(bath,1),size(bath,2)) + integer :: ilat,i,iorb,ispin_ + integer :: Nsites + logical :: check_dim + character(len=5) :: tmp_suffix + integer :: MPI_ID=0 + integer :: MPI_SIZE=1 + logical :: MPI_MASTER=.true. + ! + write(Logfile,"(A)")"" + ! + ispin_=1;if(present(ispin))ispin_=ispin + if(ispin_>Nspin)stop "ed_fit_bath_sites error: required spin index > Nspin" + ! +#ifdef _MPI + if(check_MPI())then + MPI_ID = get_Rank_MPI() + MPI_SIZE = get_Size_MPI() + MPI_MASTER = get_Master_MPI() + endif +#endif + ! + ! Check dimensions ! + Nsites=size(bath,1) + ! + do ilat = 1 + MPI_ID, Nsites, MPI_SIZE + check_dim = check_bath_dimension(bath(ilat,:)) + if(.not.check_dim) stop "init_lattice_bath: wrong bath size dimension 1 or 2 " + end do + ! call assert_shape(g,[Nsites,Nspin*Norb,Nspin*Norb,Lmats],'chi2_fitgf_generic_superc','g') do ilat=1,Nsites fg(1,ilat,:,:,:,:,:) = so2nn_reshape(g(ilat,1:Nspin*Norb,1:Nspin*Norb,1:Lmats),Nspin,Norb,Lmats) enddo - ! - rank (6) - call assert_shape(g,[Nsites,Nspin,Nspin,Norb,Norb,Lmats],'chi2_fitgf_generic_superc','g') - fg(1,:,:,:,:,:,:) = g(1:Nsites,1:Nspin,1:Nspin,1:Norb,1:Norb,1:Lmats) - end select - ! - select rank(f) - rank default - stop "chi2_fitgf_generic_superc error: input F has a wrong rank: accepted 3 or 5" - rank (3) - call assert_shape(f,[Nsites*Nspin*Norb,Nsites*Nspin*Norb,Lmats],'chi2_fitgf_generic_superc','f') - fg(2,:,:,:,:,:,:) = lso2nnn_reshape(f(1:Nsites*Nspin*Norb,1:Nsites*Nspin*Norb,1:Lmats),Nsites,Nspin,Norb,Lmats) - ! - rank (4) call assert_shape(f,[Nsites,Nspin*Norb,Nspin*Norb,Lmats],'chi2_fitgf_generic_superc','f') do ilat=1,Nsites fg(2,ilat,:,:,:,:,:) = so2nn_reshape(f(ilat,1:Nspin*Norb,1:Nspin*Norb,1:Lmats),Nspin,Norb,Lmats) enddo ! - rank (6) + bath_tmp=0.d0 + ! + do ilat= 1 + MPI_ID, Nsites, MPI_SIZE + write(Logfile,"(A)")"ed_fit_bath_sites_superc: Start Chi**2 fit for site:"//str(ilat) + bath_tmp(ilat,:) = bath(ilat,:) + call ed_set_Hloc(Hloc_ineq(ilat,:,:,:,:)) + ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) + ! + if(present(ispin))then + call ed_chi2_fitgf(fg(1,ilat,:,:,:,:,:),fg(2,ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin_,fmpi=.false.) + else + do ispin_=1,Nspin + call ed_chi2_fitgf(fg(1,ilat,:,:,:,:,:),fg(2,ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin_,fmpi=.false.) + enddo + endif + end do + ! +#ifdef _MPI + if(check_MPI())then + bath=0d0 + call AllReduce_MPI(MPI_COMM_WORLD,bath_tmp,bath) + else + bath = bath_tmp + endif +#else + bath = bath_tmp +#endif + ! + call ed_reset_suffix + end subroutine chi2_fitgf_lattice_superc_n4 + + + + + + + + subroutine chi2_fitgf_lattice_superc_n6(g,f,bath,ispin) + real(8),intent(inout) :: bath(:,:) + complex(8),dimension(:,:,:,:,:,:) :: g,f + integer,optional :: ispin + ! + complex(8) :: fg(2,size(bath,1),Nspin,Nspin,Norb,Norb,Lmats) + real(8) :: bath_tmp(size(bath,1),size(bath,2)) + integer :: ilat,i,iorb,ispin_ + integer :: Nsites + logical :: check_dim + character(len=5) :: tmp_suffix + integer :: MPI_ID=0 + integer :: MPI_SIZE=1 + logical :: MPI_MASTER=.true. + ! + write(Logfile,"(A)")"" + ! + ispin_=1;if(present(ispin))ispin_=ispin + if(ispin_>Nspin)stop "ed_fit_bath_sites error: required spin index > Nspin" + ! +#ifdef _MPI + if(check_MPI())then + MPI_ID = get_Rank_MPI() + MPI_SIZE = get_Size_MPI() + MPI_MASTER = get_Master_MPI() + endif +#endif + ! + ! Check dimensions ! + Nsites=size(bath,1) + ! + do ilat = 1 + MPI_ID, Nsites, MPI_SIZE + check_dim = check_bath_dimension(bath(ilat,:)) + if(.not.check_dim) stop "init_lattice_bath: wrong bath size dimension 1 or 2 " + end do + ! + call assert_shape(g,[Nsites,Nspin,Nspin,Norb,Norb,Lmats],'chi2_fitgf_generic_superc','g') + fg(1,:,:,:,:,:,:) = g(1:Nsites,1:Nspin,1:Nspin,1:Norb,1:Norb,1:Lmats) + ! call assert_shape(f,[Nsites,Nspin,Nspin,Norb,Norb,Lmats],'chi2_fitgf_generic_superc','f') fg(2,:,:,:,:,:,:) = f(1:Nsites,1:Nspin,1:Nspin,1:Norb,1:Norb,1:Lmats) - end select ! bath_tmp=0.d0 ! @@ -448,10 +896,10 @@ subroutine chi2_fitgf_lattice_superc(g,f,bath,ispin) ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) ! if(present(ispin))then - call ed_chi2_fitgf(fg(:,ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin_,fmpi=.false.) + call ed_chi2_fitgf(fg(1,ilat,:,:,:,:,:),fg(2,ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin_,fmpi=.false.) else do ispin_=1,Nspin - call ed_chi2_fitgf(fg(:,ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin_,fmpi=.false.) + call ed_chi2_fitgf(fg(1,ilat,:,:,:,:,:),fg(2,ilat,:,:,:,:,:),bath_tmp(ilat,:),ispin_,fmpi=.false.) enddo endif end do @@ -468,8 +916,7 @@ subroutine chi2_fitgf_lattice_superc(g,f,bath,ispin) #endif ! call ed_reset_suffix - end subroutine chi2_fitgf_lattice_superc - + end subroutine chi2_fitgf_lattice_superc_n6 diff --git a/src/ED_IO/ED_IO.f90 b/src/ED_IO/ED_IO.f90 index a0f209aa..e1fc22e7 100644 --- a/src/ED_IO/ED_IO.f90 +++ b/src/ED_IO/ED_IO.f90 @@ -13,38 +13,147 @@ MODULE ED_IO implicit none private - ! !Retrieve self-energy through routines: + ! !Retrieve self-energy through routines: interface ed_get_sigma - module procedure ed_get_sigma_site - module procedure ed_get_sigma_lattice + module procedure :: ed_get_sigma_site_n3 + module procedure :: ed_get_sigma_site_n5 + module procedure :: ed_get_sigma_lattice_n3 + module procedure :: ed_get_sigma_lattice_n4 + module procedure :: ed_get_sigma_lattice_n6 end interface ed_get_sigma !Retrieve imp GF through routines. interface ed_get_gimp - module procedure ed_get_gimp_site - module procedure ed_get_gimp_lattice + module procedure :: ed_get_gimp_site_n3 + module procedure :: ed_get_gimp_site_n5 + module procedure :: ed_get_gimp_lattice_n3 + module procedure :: ed_get_gimp_lattice_n4 + module procedure :: ed_get_gimp_lattice_n6 end interface ed_get_gimp !Retrieve imp GF_0 (G0_and) through routines. interface ed_get_g0imp - module procedure ed_get_g0imp_site - module procedure ed_get_g0imp_lattice + module procedure :: ed_get_g0imp_site_n3 + module procedure :: ed_get_g0imp_site_n5 + module procedure :: ed_get_g0imp_lattice_n3 + module procedure :: ed_get_g0imp_lattice_n4 + module procedure :: ed_get_g0imp_lattice_n6 end interface ed_get_g0imp !Rebuild impurity Sigma or GF from saved poles&weights interface ed_build_gimp - module procedure rebuild_gimp_single - module procedure rebuild_gimp_ineq + module procedure :: rebuild_gimp_single_n3 + module procedure :: rebuild_gimp_single_n5 + module procedure :: rebuild_gimp_ineq_n3 + module procedure :: rebuild_gimp_ineq_n4 + module procedure :: rebuild_gimp_ineq_n6 end interface ed_build_gimp interface ed_build_sigma - module procedure rebuild_sigma_single - module procedure rebuild_sigma_ineq + module procedure :: rebuild_sigma_single_n3 + module procedure :: rebuild_sigma_single_n5 + module procedure :: rebuild_sigma_ineq_n3 + module procedure :: rebuild_sigma_ineq_n4 + module procedure :: rebuild_sigma_ineq_n6 end interface ed_build_sigma + !Build Gand/Delta from a user bath + interface ed_get_g0and + module procedure :: ed_get_g0and_n3 + module procedure :: ed_get_g0and_n5 + end interface ed_get_g0and + interface ed_get_delta + module procedure :: ed_get_delta_n3 + module procedure :: ed_get_delta_n5 + end interface ed_get_delta + + + !Observables + interface ed_get_dens + module procedure :: ed_get_dens_n0 + module procedure :: ed_get_dens_n1 + module procedure :: ed_get_dens_n2 + end interface ed_get_dens + + interface ed_get_mag + module procedure :: ed_get_mag_n0 + module procedure :: ed_get_mag_n1 + module procedure :: ed_get_mag_n2 + module procedure :: ed_get_mag_n3 + end interface ed_get_mag + + interface ed_get_docc + module procedure :: ed_get_docc_n0 + module procedure :: ed_get_docc_n1 + module procedure :: ed_get_docc_n2 + end interface ed_get_docc + + interface ed_get_phi + module procedure :: ed_get_phisc_n0 + module procedure :: ed_get_phisc_n1 + module procedure :: ed_get_phisc_n2 + end interface ed_get_phi + + + !Get Energies + interface ed_get_eimp + module procedure :: ed_get_eimp_n1 + module procedure :: ed_get_eimp_n2 + end interface ed_get_eimp + + interface ed_get_epot + module procedure :: ed_get_epot_n0 + module procedure :: ed_get_epot_n1 + end interface ed_get_epot + + interface ed_get_eint + module procedure :: ed_get_eint_n0 + module procedure :: ed_get_eint_n1 + end interface ed_get_eint + + interface ed_get_ehartree + module procedure :: ed_get_ehartree_n0 + module procedure :: ed_get_ehartree_n1 + end interface ed_get_ehartree + + interface ed_get_eknot + module procedure :: ed_get_eknot_n0 + module procedure :: ed_get_eknot_n1 + end interface ed_get_eknot + + + !Get Double occupancies + interface ed_get_doubles + module procedure :: ed_get_doubles_n1 + module procedure :: ed_get_doubles_n2 + end interface ed_get_doubles + + interface ed_get_dust + module procedure :: ed_get_dust_n0 + module procedure :: ed_get_dust_n1 + end interface ed_get_dust + + interface ed_get_dund + module procedure :: ed_get_dund_n0 + module procedure :: ed_get_dund_n1 + end interface ed_get_dund + + interface ed_get_dse + module procedure :: ed_get_dse_n0 + module procedure :: ed_get_dse_n1 + end interface ed_get_dse + + interface ed_get_dph + module procedure :: ed_get_dph_n0 + module procedure :: ed_get_dph_n1 + end interface ed_get_dph + + + + interface ed_get_density_matrix module procedure :: ed_get_density_matrix_single @@ -69,7 +178,7 @@ MODULE ED_IO public :: ed_get_dens public :: ed_get_mag public :: ed_get_docc - public :: ed_get_phisc + public :: ed_get_phi public :: ed_get_eimp public :: ed_get_epot public :: ed_get_eint @@ -153,7 +262,13 @@ MODULE ED_IO !+--------------------------------------------------------------------------+! ! PURPOSE: Retrieve measured values of the local observables !+--------------------------------------------------------------------------+! - include "get_observables.f90" + include "get_dens.f90" + include "get_mag.f90" + include "get_docc.f90" + include "get_phi.f90" + include "get_energy.f90" + include "get_doubles.f90" + include "get_neigen.f90" diff --git a/src/ED_IO/get_dens.f90 b/src/ED_IO/get_dens.f90 new file mode 100644 index 00000000..0896076d --- /dev/null +++ b/src/ED_IO/get_dens.f90 @@ -0,0 +1,38 @@ +subroutine ed_get_dens_n0(self,iorb) + real(8) :: self + integer,optional :: iorb + integer :: iorb_ + iorb_=1;if(present(iorb))iorb_=iorb + if(iorb_>Norb)stop "ed_get_dens error: orbital index > N_orbital" + self = ed_dens(iorb_) +end subroutine ed_get_dens_n0 + + +subroutine ed_get_dens_n1(self,iorb,Nlat) + real(8),dimension(:) :: self + integer,optional :: iorb,Nlat + integer :: iorb_ + iorb_=1;if(present(iorb))iorb_=iorb + if(iorb_>Norb)stop "ed_get_dens error: orbital index > N_orbital" + if(present(Nlat))then + if(.not.allocated(dens_ineq))stop "ed_get_dens error: dens_ineq not allocated" + if(Nlat>size(dens_ineq,1))stop "ed_get_dens error: required N_sites > evaluated N_sites" + endif + if(present(Nlat))then + call assert_shape(self,[Nlat],'ed_get_dens','dens') + self = dens_ineq(:,iorb_) + else + call assert_shape(self,[Norb],'ed_get_dens','dens') + self = ed_dens + endif +end subroutine ed_get_dens_n1 + + +subroutine ed_get_dens_n2(self,Nlat) + real(8),dimension(:,:) :: self + integer :: Nlat + if(.not.allocated(dens_ineq))stop "ed_get_dens error: dens_ineq not allocated" + if(Nlat>size(dens_ineq,1))stop "ed_get_dens error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat,Norb],'ed_get_dens','dens') + self = dens_ineq +end subroutine ed_get_dens_n2 diff --git a/src/ED_IO/get_docc.f90 b/src/ED_IO/get_docc.f90 new file mode 100644 index 00000000..53aea7ea --- /dev/null +++ b/src/ED_IO/get_docc.f90 @@ -0,0 +1,36 @@ +subroutine ed_get_docc_n0(self,iorb) + real(8) :: self + integer,optional :: iorb + integer :: iorb_ + iorb_=1;if(present(iorb))iorb_=iorb + if(iorb_>Norb)stop "ed_get_docc error: orbital index > N_orbital" + self = ed_docc(iorb_) +end subroutine ed_get_docc_n0 + +subroutine ed_get_docc_n1(self,iorb,Nlat) + real(8),dimension(:) :: self + integer,optional :: iorb,Nlat + integer :: iorb_ + iorb_=1;if(present(iorb))iorb_=iorb + if(iorb_>Norb)stop "ed_get_docc error: orbital index > N_orbital" + if(present(Nlat))then + if(.not.allocated(docc_ineq))stop "ed_get_docc error: docc_ineq not allocated" + if(Nlat>size(docc_ineq,1))stop "ed_get_docc error: required N_sites > evaluated N_sites" + endif + if(present(Nlat))then + call assert_shape(self,[Nlat],'ed_get_docc','docc') + self = docc_ineq(:,iorb_) + else + call assert_shape(self,[Norb],'ed_get_docc','docc') + self = ed_docc + endif +end subroutine ed_get_docc_n1 + +subroutine ed_get_docc_n2(self,Nlat) + real(8),dimension(:,:) :: self + integer :: Nlat + if(.not.allocated(docc_ineq))stop "ed_get_docc error: docc_ineq not allocated" + if(Nlat>size(docc_ineq,1))stop "ed_get_docc error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat,Norb],'ed_get_docc','docc') + self = docc_ineq +end subroutine ed_get_docc_n2 diff --git a/src/ED_IO/get_doubles.f90 b/src/ED_IO/get_doubles.f90 new file mode 100644 index 00000000..1e70b4cf --- /dev/null +++ b/src/ED_IO/get_doubles.f90 @@ -0,0 +1,88 @@ +subroutine ed_get_doubles_n1(self) + real(8),dimension(:) :: self + call assert_shape(self,[4],'ed_get_doubles','doubles') + self = [ed_Dust,ed_Dund,ed_Dse,ed_Dph] +end subroutine ed_get_doubles_n1 + +subroutine ed_get_doubles_n2(self,Nlat) + real(8),dimension(:,:) :: self + integer :: Nlat + ! + if(.not.allocated(dd_ineq))stop "ed_get_doubles error: dd_ineq not allocated" + if(Nlat>size(dd_ineq,1))stop "ed_get_doubles error: required N_sites > evaluated N_sites" + call assert_shape(self,[2,Nlat],'ed_get_doubles','doubles') + self = dd_ineq +end subroutine ed_get_doubles_n2 + + + + + +subroutine ed_get_dust_n0(self) + real(8) :: self + self = ed_Dust +end subroutine ed_get_dust_n0 + +subroutine ed_get_dust_n1(self,Nlat) + real(8),dimension(:) :: self + integer :: Nlat + ! + if(.not.allocated(dd_ineq))stop "ed_get_dust error: dd_ineq not allocated" + if(Nlat>size(dd_ineq,1))stop "ed_get_dust error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat],'ed_get_dust','dust') + self = dd_ineq(:,1) +end subroutine ed_get_dust_n1 + + + + + +subroutine ed_get_dund_n0(self) + real(8) :: self + self = ed_Dund +end subroutine ed_get_dund_n0 + +subroutine ed_get_dund_n1(self,Nlat) + real(8),dimension(:) :: self + integer :: Nlat + ! + if(.not.allocated(dd_ineq))stop "ed_get_dund error: dd_ineq not allocated" + if(Nlat>size(dd_ineq,1))stop "ed_get_dund error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat],'ed_get_dund','dund') + self = dd_ineq(:,2) +end subroutine ed_get_dund_n1 + + + +subroutine ed_get_dse_n0(self) + real(8) :: self + self = ed_Dse +end subroutine ed_get_dse_n0 + +subroutine ed_get_dse_n1(self,Nlat) + real(8),dimension(:) :: self + integer :: Nlat + ! + if(.not.allocated(dd_ineq))stop "ed_get_dse error: dd_ineq not allocated" + if(Nlat>size(dd_ineq,1))stop "ed_get_dse error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat],'ed_get_dse','dse') + self = dd_ineq(:,3) +end subroutine ed_get_dse_n1 + + + + +subroutine ed_get_dph_n0(self) + real(8) :: self + self = ed_Dph +end subroutine ed_get_dph_n0 + +subroutine ed_get_dph_n1(self,Nlat) + real(8),dimension(:) :: self + integer :: Nlat + ! + if(.not.allocated(dd_ineq))stop "ed_get_dph error: dd_ineq not allocated" + if(Nlat>size(dd_ineq,1))stop "ed_get_dph error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat],'ed_get_dph','dph') + self = dd_ineq(:,4) +end subroutine ed_get_dph_n1 diff --git a/src/ED_IO/get_energy.f90 b/src/ED_IO/get_energy.f90 new file mode 100644 index 00000000..26cb9c13 --- /dev/null +++ b/src/ED_IO/get_energy.f90 @@ -0,0 +1,90 @@ +subroutine ed_get_eimp_n1(self) + real(8),dimension(:) :: self + call assert_shape(self,[4],'ed_get_eimp','eimp') + self = [ed_Epot,ed_Eint,ed_Ehartree,ed_Eknot] +end subroutine ed_get_eimp_n1 + +subroutine ed_get_eimp_n2(self,Nlat) + real(8),dimension(:,:) :: self + integer :: Nlat + ! + if(.not.allocated(e_ineq))stop "ed_get_eimp error: e_ineq not allocated" + if(Nlat>size(e_ineq,1))stop "ed_get_eimp error: required N_sites > evaluated N_sites" + call assert_shape(self,[2,Nlat],'ed_get_eimp','eimp') + self = e_ineq +end subroutine ed_get_eimp_n2 + + + + +subroutine ed_get_epot_n0(self) + real(8) :: self + self = ed_Epot +end subroutine ed_get_epot_n0 + +subroutine ed_get_epot_n1(self,Nlat) + real(8),dimension(:) :: self + integer :: Nlat + ! + if(.not.allocated(e_ineq))stop "ed_get_epot error: e_ineq not allocated" + if(Nlat>size(e_ineq,1))stop "ed_get_epot error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat],'ed_get_epot','epot') + self = e_ineq(:,1) +end subroutine ed_get_epot_n1 + + + + +subroutine ed_get_eint_n0(self) + real(8) :: self + self = ed_Eint +end subroutine ed_get_eint_n0 + +subroutine ed_get_eint_n1(self,Nlat) + real(8),dimension(:) :: self + integer :: Nlat + ! + if(.not.allocated(e_ineq))stop "ed_get_eint error: e_ineq not allocated" + if(Nlat>size(e_ineq,1))stop "ed_get_eint error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat],'ed_get_eint','eint') + self = e_ineq(:,2) +end subroutine ed_get_eint_n1 + + + + + +subroutine ed_get_ehartree_n0(self) + real(8) :: self + self = ed_Ehartree +end subroutine ed_get_ehartree_n0 + +subroutine ed_get_ehartree_n1(self,Nlat) + real(8),dimension(:) :: self + integer :: Nlat + ! + if(.not.allocated(e_ineq))stop "ed_get_ehartree error: e_ineq not allocated" + if(Nlat>size(e_ineq,1))stop "ed_get_ehartree error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat],'ed_get_ehartree','ehartree') + self = e_ineq(:,3) +end subroutine ed_get_ehartree_n1 + + + + +subroutine ed_get_eknot_n0(self) + real(8) :: self + self = ed_Eknot +end subroutine ed_get_eknot_n0 + +subroutine ed_get_eknot_n1(self,Nlat) + real(8),dimension(:) :: self + integer :: Nlat + ! + if(.not.allocated(e_ineq))stop "ed_get_eknot error: e_ineq not allocated" + if(Nlat>size(e_ineq,1))stop "ed_get_eknot error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat],'ed_get_eknot','eknot') + self = e_ineq(:,4) +end subroutine ed_get_eknot_n1 + + diff --git a/src/ED_IO/get_g0imp.f90 b/src/ED_IO/get_g0imp.f90 index 6f90b6ee..0e8295de 100644 --- a/src/ED_IO/get_g0imp.f90 +++ b/src/ED_IO/get_g0imp.f90 @@ -1,5 +1,5 @@ -subroutine ed_get_g0imp_site(self,bath,axis,type) - complex(8),dimension(..),intent(inout) :: self +subroutine ed_get_g0imp_site_n3(self,bath,axis,type) + complex(8),dimension(:,:,:),intent(inout) :: self real(8),dimension(:) :: bath character(len=*),optional :: axis character(len=*),optional :: type @@ -27,21 +27,53 @@ subroutine ed_get_g0imp_site(self,bath,axis,type) allocate(g0(Nspin,Nspin,Norb,Norb,L)) call ed_get_g0and(x,bath,g0,axis=axis_,type=type_) ! - select rank(self) - rank default; stop "ed_get_g0imp ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,L],'ed_get_g0imp','self') self = nn2so_reshape(g0,Nspin,Norb,L) - rank (5) +end subroutine ed_get_g0imp_site_n3 + +subroutine ed_get_g0imp_site_n5(self,bath,axis,type) + complex(8),dimension(:,:,:,:,:),intent(inout) :: self + real(8),dimension(:) :: bath + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + complex(8),dimension(:,:,:,:,:),allocatable :: g0 + complex(8),dimension(:),allocatable :: x + integer :: L + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + call allocate_grids + ! + select case (axis_) + case default;stop "ed_get_g0imp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + L = Lmats + allocate(x(L)) + x=xi*wm + case('r','R') + L = Lreal + allocate(x(L)) + x=dcmplx(wr,eps) + end select + ! + allocate(g0(Nspin,Nspin,Norb,Norb,L)) + call ed_get_g0and(x,bath,g0,axis=axis_,type=type_) + ! call assert_shape(self,[Nspin,Nspin,Norb,Norb,L],'ed_get_g0imp','self') self = g0 - end select -end subroutine ed_get_g0imp_site +end subroutine ed_get_g0imp_site_n5 + + + +!################################################################## +!################################################################## +!################################################################## -subroutine ed_get_g0imp_lattice(self,bath,axis,type) - complex(8),dimension(..),intent(inout) :: self +subroutine ed_get_g0imp_lattice_n3(self,bath,axis,type) + complex(8),dimension(:,:,:),intent(inout) :: self real(8),dimension(:,:) :: bath character(len=*),optional :: axis character(len=*),optional :: type @@ -53,7 +85,7 @@ subroutine ed_get_g0imp_lattice(self,bath,axis,type) axis_='m';if(present(axis))axis_=trim(axis) type_='n';if(present(type))type_=trim(type) call allocate_grids - select case (axis_) + select case (axis_) case default;stop "ed_get_g0imp ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') L = Lmats @@ -71,18 +103,79 @@ subroutine ed_get_g0imp_lattice(self,bath,axis,type) call ed_get_g0and(x,bath(isite,:),g0(isite,:,:,:,:,:),axis=axis_,type=type_) enddo ! - select rank(self) - rank default; stop "ed_get_g0imp ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nsites*Nspin*Norb,Nsites*Nspin*Norb,L],'ed_get_g0imp','self') self = nnn2lso_reshape(g0,Nsites,Nspin,Norb,L) - rank (4) +end subroutine ed_get_g0imp_lattice_n3 + +subroutine ed_get_g0imp_lattice_n4(self,bath,axis,type) + complex(8),dimension(:,:,:,:),intent(inout) :: self + real(8),dimension(:,:) :: bath + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + complex(8),dimension(:,:,:,:,:,:),allocatable :: g0 + complex(8),dimension(:),allocatable :: x + integer :: L,Nsites,isite + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + call allocate_grids + select case (axis_) + case default;stop "ed_get_g0imp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + L = Lmats + allocate(x(L)) + x=xi*wm + case('r','R') + L = Lreal + allocate(x(L)) + x=dcmplx(wr,eps) + end select + ! + Nsites=size(bath,1) + allocate(g0(Nsites,Nspin,Nspin,Norb,Norb,L)) + do isite=1,Nsites + call ed_get_g0and(x,bath(isite,:),g0(isite,:,:,:,:,:),axis=axis_,type=type_) + enddo + ! call assert_shape(self,[Nsites,Nspin*Norb,Nspin*Norb,L],'ed_get_g0imp','self') do isite=1,Nsites self(isite,:,:,:) = nn2so_reshape(g0(isite,:,:,:,:,:),Nspin,Norb,L) enddo - rank (6) +end subroutine ed_get_g0imp_lattice_n4 + + +subroutine ed_get_g0imp_lattice_n6(self,bath,axis,type) + complex(8),dimension(:,:,:,:,:,:),intent(inout) :: self + real(8),dimension(:,:) :: bath + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + complex(8),dimension(:,:,:,:,:,:),allocatable :: g0 + complex(8),dimension(:),allocatable :: x + integer :: L,Nsites,isite + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + call allocate_grids + select case (axis_) + case default;stop "ed_get_g0imp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + L = Lmats + allocate(x(L)) + x=xi*wm + case('r','R') + L = Lreal + allocate(x(L)) + x=dcmplx(wr,eps) + end select + ! + Nsites=size(bath,1) + allocate(g0(Nsites,Nspin,Nspin,Norb,Norb,L)) + do isite=1,Nsites + call ed_get_g0and(x,bath(isite,:),g0(isite,:,:,:,:,:),axis=axis_,type=type_) + enddo + ! call assert_shape(self,[Nsites,Nspin,Nspin,Norb,Norb,L],'ed_get_g0imp','self') self = g0 - end select -end subroutine ed_get_g0imp_lattice +end subroutine ed_get_g0imp_lattice_n6 diff --git a/src/ED_IO/get_gand_bath.f90 b/src/ED_IO/get_gand_bath.f90 index affdc784..f5db3946 100644 --- a/src/ED_IO/get_gand_bath.f90 +++ b/src/ED_IO/get_gand_bath.f90 @@ -1,7 +1,7 @@ -subroutine ed_get_g0and(x,bath_,G0and,axis,type) +subroutine ed_get_g0and_n3(x,bath_,G0and,axis,type) complex(8),dimension(:),intent(in) :: x real(8),dimension(:) :: bath_ - complex(8),dimension(..) :: G0and + complex(8),dimension(:,:,:) :: G0and character(len=*),optional :: axis character(len=*),optional :: type ! @@ -27,23 +27,52 @@ subroutine ed_get_g0and(x,bath_,G0and,axis,type) call deallocate_dmft_bath(dmft_bath_) ! L=size(x) - select rank(g0and) - rank default; stop "ed_get_g0and ERROR: g0and has a wrong rank" - rank (3) call assert_shape(g0and,[Nspin*Norb,Nspin*Norb,L],'ed_get_g0and','g0and') g0and = nn2so_reshape(g0,Nspin,Norb,L) - rank (5) +end subroutine ed_get_g0and_n3 + +subroutine ed_get_g0and_n5(x,bath_,G0and,axis,type) + complex(8),dimension(:),intent(in) :: x + real(8),dimension(:) :: bath_ + complex(8),dimension(:,:,:,:,:) :: G0and + character(len=*),optional :: axis + character(len=*),optional :: type + ! + type(effective_bath) :: dmft_bath_ + logical :: check + character(len=1) :: axis_,type_ + complex(8),dimension(Nspin,Nspin,Norb,Norb,size(x)) :: g0 + integer :: L + ! + axis_='m';if(present(axis))axis_=axis + type_='n';if(present(type))type_=trim(type) + check= check_bath_dimension(bath_) + if(.not.check)stop "g0and_bath_mats_main_ error: wrong bath dimensions" + call allocate_dmft_bath(dmft_bath_) + call set_dmft_bath(bath_,dmft_bath_) + select case(type_) + case default;stop "ed_get_g0and ERROR: type is wrong: either Normal or Anomalous" + case ('n','N') + g0 = g0and_bath_function(x,dmft_bath_) + case('a','A') + g0 = f0and_bath_function(x,dmft_bath_) + end select + call deallocate_dmft_bath(dmft_bath_) + ! + L=size(x) call assert_shape(g0and,[Nspin,Nspin,Norb,Norb,L],'ed_get_g0and','g0and') g0and = g0 - end select -end subroutine ed_get_g0and +end subroutine ed_get_g0and_n5 + -subroutine ed_get_delta(x,bath_,delta,axis,type) + + +subroutine ed_get_delta_n3(x,bath_,delta,axis,type) complex(8),dimension(:),intent(in) :: x real(8),dimension(:) :: bath_ - complex(8),dimension(..) :: delta + complex(8),dimension(:,:,:) :: delta character(len=*),optional :: axis character(len=*),optional :: type ! @@ -68,16 +97,41 @@ subroutine ed_get_delta(x,bath_,delta,axis,type) call deallocate_dmft_bath(dmft_bath_) ! L=size(x) - select rank(delta) - rank default; stop "ed_get_delta ERROR: delta has a wrong rank" - rank (3) call assert_shape(delta,[Nspin*Norb,Nspin*Norb,L],'ed_get_delta','delta') delta = nn2so_reshape(d0,Nspin,Norb,L) - rank (5) +end subroutine ed_get_delta_n3 + +subroutine ed_get_delta_n5(x,bath_,delta,axis,type) + complex(8),dimension(:),intent(in) :: x + real(8),dimension(:) :: bath_ + complex(8),dimension(:,:,:,:,:) :: delta + character(len=*),optional :: axis + character(len=*),optional :: type + ! + type(effective_bath) :: dmft_bath_ + logical :: check + character(len=1) :: axis_,type_ + complex(8),dimension(Nspin,Nspin,Norb,Norb,size(x)) :: d0 + integer :: L + ! + axis_='mats';if(present(axis))axis_=axis + check= check_bath_dimension(bath_) + if(.not.check)stop "delta_bath_mats_main_ error: wrong bath dimensions" + call allocate_dmft_bath(dmft_bath_) + call set_dmft_bath(bath_,dmft_bath_) + select case(type_) + case default;stop "ed_get_delta ERROR: type is wrong: either Normal or Anomalous" + case ('n','N') + d0 = delta_bath_function(x,dmft_bath_,axis_) + case('a','A') + d0 = fdelta_bath_function(x,dmft_bath_,axis_) + end select + call deallocate_dmft_bath(dmft_bath_) + ! + L=size(x) call assert_shape(delta,[Nspin,Nspin,Norb,Norb,L],'ed_get_delta','delta') delta = d0 - end select -end subroutine ed_get_delta +end subroutine ed_get_delta_n5 diff --git a/src/ED_IO/get_gimp.f90 b/src/ED_IO/get_gimp.f90 index 295198a4..ead4b8dc 100644 --- a/src/ED_IO/get_gimp.f90 +++ b/src/ED_IO/get_gimp.f90 @@ -1,9 +1,9 @@ -subroutine ed_get_gimp_site(self,axis,type) - complex(8),dimension(..),intent(inout) :: self - character(len=*),optional :: axis - character(len=*),optional :: type - character(len=1) :: axis_ - character(len=1) :: type_ +subroutine ed_get_gimp_site_n3(self,axis,type) + complex(8),dimension(:,:,:),intent(inout) :: self + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ axis_='m';if(present(axis))axis_=trim(axis) type_='n';if(present(type))type_=trim(type) select case(type_) @@ -12,57 +12,67 @@ subroutine ed_get_gimp_site(self,axis,type) select case(axis_) case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') - select rank(self) - rank default; stop "ed_get_gimp ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,Lmats],'ed_get_gimp','self') self = nn2so_reshape(impGmats,Nspin,Norb,Lmats) - rank (5) - call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lmats],'ed_get_gimp','self') - self = impGmats - end select case('r','R') - select rank(self) - rank default; stop "ed_get_gimp ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,Lreal],'ed_get_gimp','self') self = nn2so_reshape(impGreal,Nspin,Norb,Lmats) - rank (5) - call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lreal],'ed_get_gimp','self') - self = impGreal - end select end select case('a','A') select case(axis_) case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') - select rank(self) - rank default; stop "ed_get_gimp ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,Lmats],'ed_get_gimp','self') self = nn2so_reshape(impFmats,Nspin,Norb,Lmats) - rank (5) - call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lmats],'ed_get_gimp','self') - self = impFmats - end select case('r','R') - select rank(self) - rank default; stop "ed_get_gimp ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,Lreal],'ed_get_gimp','self') self = nn2so_reshape(impFreal,Nspin,Norb,Lmats) - rank (5) + end select + end select +end subroutine ed_get_gimp_site_n3 + + +subroutine ed_get_gimp_site_n5(self,axis,type) + complex(8),dimension(:,:,:,:,:),intent(inout) :: self + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + select case(type_) + case default; stop "ed_get_gimp ERROR: type is neither Normal, nor Anomalous" + case ('n','N') + select case(axis_) + case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lmats],'ed_get_gimp','self') + self = impGmats + case('r','R') + call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lreal],'ed_get_gimp','self') + self = impGreal + end select + case('a','A') + select case(axis_) + case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lmats],'ed_get_gimp','self') + self = impFmats + case('r','R') call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lreal],'ed_get_gimp','self') self = impFreal - end select end select end select -end subroutine ed_get_gimp_site +end subroutine ed_get_gimp_site_n5 + +!################################################################## +!################################################################## +!################################################################## -subroutine ed_get_gimp_lattice(self,nlat,axis,type) - complex(8),dimension(..),intent(inout) :: self +subroutine ed_get_gimp_lattice_n3(self,nlat,axis,type) + complex(8),dimension(:,:,:),intent(inout) :: self integer,intent(in) :: nlat character(len=*),optional :: axis character(len=*),optional :: type @@ -78,72 +88,104 @@ subroutine ed_get_gimp_lattice(self,nlat,axis,type) select case(axis_) case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') - select rank(self) - rank default; stop "ed_get_gimp ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lmats],'ed_get_gimp','self') self = nnn2lso_reshape(Gmats_ineq,Nlat,Nspin,Norb,Lmats) - rank (4) + case('r','R') + call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lreal],'ed_get_gimp','self') + self = nnn2lso_reshape(Freal_ineq,Nlat,Nspin,Norb,Lreal) + end select + case('a','A') + select case(axis_) + case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lmats],'ed_get_gimp','self') + self = nnn2lso_reshape(Fmats_ineq,Nlat,Nspin,Norb,Lmats) + case('r','R') + call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lreal],'ed_get_gimp','self') + self = nnn2lso_reshape(Freal_ineq,Nlat,Nspin,Norb,Lreal) + end select + end select +end subroutine ed_get_gimp_lattice_n3 + +subroutine ed_get_gimp_lattice_n4(self,nlat,axis,type) + complex(8),dimension(:,:,:,:),intent(inout) :: self + integer,intent(in) :: nlat + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + integer :: ilat + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + if(Nlat/=size(Gmats_ineq,1))stop "ERROR ed_get_gimp: wrong Nlat" + select case(type_) + case default; stop "ed_get_gimp ERROR: type is neither Normal, nor Anomalous" + case ('n','N') + select case(axis_) + case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') call assert_shape(self,[Nlat,Nspin*Norb,Nspin*Norb,Lmats],'ed_get_gimp','self') do ilat=1,Nlat self(ilat,:,:,:) = nn2so_reshape(Gmats_ineq(ilat,:,:,:,:,:),Nspin,Norb,Lmats) enddo - rank (6) - call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lmats],'ed_get_gimp','self') - self = Gmats_ineq - end select case('r','R') - select rank(self) - rank default; stop "ed_get_gimp ERROR: self has a wrong rank" - rank (3) - call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lreal],'ed_get_gimp','self') - self = nnn2lso_reshape(Freal_ineq,Nlat,Nspin,Norb,Lreal) - rank (4) call assert_shape(self,[Nlat,Nspin*Norb,Nspin*Norb,Lreal],'ed_get_gimp','self') do ilat=1,Nlat self(ilat,:,:,:) = nn2so_reshape(Freal_ineq(ilat,:,:,:,:,:),Nspin,Norb,Lreal) enddo - rank (6) - call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lreal],'ed_get_gimp','self') - self = Freal_ineq - end select end select case('a','A') select case(axis_) case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') - select rank(self) - rank default; stop "ed_get_gimp ERROR: self has a wrong rank" - rank (3) - call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lmats],'ed_get_gimp','self') - self = nnn2lso_reshape(Fmats_ineq,Nlat,Nspin,Norb,Lmats) - rank (4) call assert_shape(self,[Nlat,Nspin*Norb,Nspin*Norb,Lmats],'ed_get_gimp','self') do ilat=1,Nlat self(ilat,:,:,:) = nn2so_reshape(Fmats_ineq(ilat,:,:,:,:,:),Nspin,Norb,Lmats) enddo - rank (6) - call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lmats],'ed_get_gimp','self') - self = Fmats_ineq - end select case('r','R') - select rank(self) - rank default; stop "ed_get_gimp ERROR: self has a wrong rank" - rank (3) - call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lreal],'ed_get_gimp','self') - self = nnn2lso_reshape(Freal_ineq,Nlat,Nspin,Norb,Lreal) - rank (4) call assert_shape(self,[Nlat,Nspin*Norb,Nspin*Norb,Lreal],'ed_get_gimp','self') do ilat=1,Nlat self(ilat,:,:,:) = nn2so_reshape(Freal_ineq(ilat,:,:,:,:,:),Nspin,Norb,Lreal) enddo - rank (6) + end select + end select +end subroutine ed_get_gimp_lattice_n4 + +subroutine ed_get_gimp_lattice_n6(self,nlat,axis,type) + complex(8),dimension(:,:,:,:,:,:),intent(inout) :: self + integer,intent(in) :: nlat + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + integer :: ilat + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + if(Nlat/=size(Gmats_ineq,1))stop "ERROR ed_get_gimp: wrong Nlat" + select case(type_) + case default; stop "ed_get_gimp ERROR: type is neither Normal, nor Anomalous" + case ('n','N') + select case(axis_) + case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lmats],'ed_get_gimp','self') + self = Gmats_ineq + case('r','R') + call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lreal],'ed_get_gimp','self') + self = Freal_ineq + end select + case('a','A') + select case(axis_) + case default;stop "ed_get_gimp ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lmats],'ed_get_gimp','self') + self = Fmats_ineq + case('r','R') call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lreal],'ed_get_gimp','self') self = Freal_ineq - end select end select end select -end subroutine ed_get_gimp_lattice +end subroutine ed_get_gimp_lattice_n6 diff --git a/src/ED_IO/get_mag.f90 b/src/ED_IO/get_mag.f90 new file mode 100644 index 00000000..1877a027 --- /dev/null +++ b/src/ED_IO/get_mag.f90 @@ -0,0 +1,89 @@ +subroutine ed_get_mag_n0(self,component,iorb) + real(8) :: self + character(len=1),optional :: component + integer,optional :: iorb + ! + integer :: iorb_ + character(len=1) :: char + integer :: id + ! + iorb_=1;if(present(iorb))iorb_=iorb + char='Z';if(present(component))char=component + select case(char) + case default;id=3 + case('X','x');id=1 + case('Y','y');id=2 + end select + if(iorb_>Norb)stop "ed_get_mag error: orbital index > N_orbital" + self = ed_mag(id,iorb_) +end subroutine ed_get_mag_n0 + +subroutine ed_get_mag_n1(self,component,iorb,Nlat) + real(8),dimension(:) :: self + character(len=1),optional :: component + integer,optional :: iorb,Nlat + ! + integer :: iorb_ + character(len=1) :: char + integer :: id + ! + iorb_=1;if(present(iorb))iorb_=iorb + char='Z';if(present(component))char=component + select case(char) + case default;id=3 + case('X','x');id=1 + case('Y','y');id=2 + end select + if(iorb_>Norb)stop "ed_get_mag error: orbital index > N_orbital" + if(present(Nlat))then + if(.not.allocated(mag_ineq))stop "ed_get_mag error: mag_ineq not allocated" + if(Nlat>size(mag_ineq,1))stop "ed_get_mag error: required N_sites > evaluated N_sites" + endif + if(present(Nlat))then + call assert_shape(self,[Nlat],'ed_get_mag','mag') + self = mag_ineq(:,id,iorb_) + else + call assert_shape(self,[Norb],'ed_get_mag','mag') + self = ed_mag(id,:) + endif +end subroutine ed_get_mag_n1 + +subroutine ed_get_mag_n2(self,component,Nlat) + real(8),dimension(:,:) :: self + character(len=1),optional :: component + integer :: Nlat + ! + character(len=1) :: char + integer :: id + ! + char='Z';if(present(component))char=component + select case(char) + case default;id=3 + case('X','x');id=1 + case('Y','y');id=2 + end select + if(.not.allocated(mag_ineq))stop "ed_get_mag error: mag_ineq not allocated" + if(Nlat>size(mag_ineq,1))stop "ed_get_mag error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat,Norb],'ed_get_mag','mag') + self = mag_ineq(:,id,:) +end subroutine ed_get_mag_n2 + +subroutine ed_get_mag_n3(self,component,Nlat) + real(8),dimension(:,:,:) :: self + character(len=1),optional :: component + integer :: Nlat + ! + character(len=1) :: char + integer :: id + ! + char='Z';if(present(component))char=component + select case(char) + case default;id=3 + case('X','x');id=1 + case('Y','y');id=2 + end select + if(.not.allocated(mag_ineq))stop "ed_get_mag error: mag_ineq not allocated" + if(Nlat>size(mag_ineq,1))stop "ed_get_mag error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat,3,Norb],'ed_get_mag','mag') + self = mag_ineq +end subroutine ed_get_mag_n3 diff --git a/src/ED_IO/get_neigen.f90 b/src/ED_IO/get_neigen.f90 new file mode 100644 index 00000000..14076d99 --- /dev/null +++ b/src/ED_IO/get_neigen.f90 @@ -0,0 +1,11 @@ +subroutine ed_get_neigen_total(nlii,Nlat) + integer :: Nlat + integer,dimension(Nlat) :: nlii + nlii=0d0 + if(allocated(neigen_total_ineq))then + if(Nlat>size(neigen_total_ineq)) stop "ed_get_docc error: required N_sites > evaluated N_sites" + nlii=neigen_total_ineq + endif +end subroutine ed_get_neigen_total + + diff --git a/src/ED_IO/get_observables.f90 b/src/ED_IO/get_observables.f90 deleted file mode 100644 index 0a0bc8dd..00000000 --- a/src/ED_IO/get_observables.f90 +++ /dev/null @@ -1,359 +0,0 @@ -subroutine ed_get_dens(self,iorb,Nlat) - real(8),dimension(..) :: self - integer,optional :: iorb,Nlat - integer :: iorb_ - iorb_=1;if(present(iorb))iorb_=iorb - if(iorb_>Norb)stop "ed_get_dens error: orbital index > N_orbital" - if(present(Nlat))then - if(.not.allocated(dens_ineq))stop "ed_get_dens error: dens_ineq not allocated" - if(Nlat>size(dens_ineq,1))stop "ed_get_dens error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_dens error: dens has wrong rank' - rank (0) !scalar - self = ed_dens(iorb_) - rank (1) !all orbitals or all sites 1orbital - if(present(Nlat))then - call assert_shape(self,[Nlat],'ed_get_dens','dens') - self = dens_ineq(:,iorb_) - else - call assert_shape(self,[Norb],'ed_get_dens','dens') - self = ed_dens - endif - rank (2) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_dens error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat,Norb],'ed_get_dens','dens') - self = dens_ineq - end select -end subroutine ed_get_dens - - - -subroutine ed_get_mag(self,component,iorb,Nlat) - real(8),dimension(..) :: self - character(len=1),optional :: component - integer,optional :: iorb,Nlat - ! - integer :: iorb_ - character(len=1) :: char - integer :: id - ! - iorb_=1;if(present(iorb))iorb_=iorb - char='Z';if(present(component))char=component - select case(char) - case default;id=3 - case('X','x');id=1 - case('Y','y');id=2 - end select - if(iorb_>Norb)stop "ed_get_mag error: orbital index > N_orbital" - if(present(Nlat))then - if(.not.allocated(mag_ineq))stop "ed_get_mag error: mag_ineq not allocated" - if(Nlat>size(mag_ineq,1))stop "ed_get_mag error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_mag error: mag has wrong rank' - rank (0) !scalar - self = ed_mag(id,iorb_) - rank (1) - if(present(Nlat))then - call assert_shape(self,[Nlat],'ed_get_mag','mag') - self = mag_ineq(:,id,iorb_) - else - call assert_shape(self,[Norb],'ed_get_mag','mag') - self = ed_mag(id,:) - endif - rank (2) - if(.not.present(Nlat))stop "ed_get_mag error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat,Norb],'ed_get_mag','mag') - self = mag_ineq(:,id,:) - rank (3) - if(.not.present(Nlat))stop "ed_get_mag error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat,3,Norb],'ed_get_mag','mag') - self = mag_ineq - end select -end subroutine ed_get_mag - - -subroutine ed_get_docc(self,iorb,Nlat) - real(8),dimension(..) :: self - integer,optional :: iorb,Nlat - integer :: iorb_ - iorb_=1;if(present(iorb))iorb_=iorb - if(iorb_>Norb)stop "ed_get_docc error: orbital index > N_orbital" - if(present(Nlat))then - if(.not.allocated(docc_ineq))stop "ed_get_docc error: docc_ineq not allocated" - if(Nlat>size(docc_ineq,1))stop "ed_get_docc error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_docc error: docc has wrong rank' - rank (0) !scalar - self = ed_docc(iorb_) - rank (1) !all orbitals or all sites 1orbital - if(present(Nlat))then - call assert_shape(self,[Nlat],'ed_get_docc','docc') - self = docc_ineq(:,iorb_) - else - call assert_shape(self,[Norb],'ed_get_docc','docc') - self = ed_docc - endif - rank (2) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_docc error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat,Norb],'ed_get_docc','docc') - self = docc_ineq - end select -end subroutine ed_get_docc - - -subroutine ed_get_phisc(self,iorb,Nlat) - real(8),dimension(..) :: self - integer,optional :: iorb,Nlat - integer :: iorb_ - iorb_=1;if(present(iorb))iorb_=iorb - if(iorb_>Norb)stop "ed_get_phisc error: orbital index > N_orbital" - if(present(Nlat))then - if(.not.allocated(phisc_ineq))stop "ed_get_phisc error: phisc_ineq not allocated" - if(Nlat>size(phisc_ineq,1))stop "ed_get_phisc error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_phisc error: phisc has wrong rank' - rank (0) !scalar - self = ed_phisc(iorb_) - rank (1) !all orbitals or all sites 1orbital - if(present(Nlat))then - call assert_shape(self,[Nlat],'ed_get_phisc','phisc') - self = phisc_ineq(:,iorb_) - else - call assert_shape(self,[Norb],'ed_get_phisc','phisc') - self = ed_phisc - endif - rank (2) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_phisc error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat,Norb],'ed_get_phisc','phisc') - self = phisc_ineq - end select -end subroutine ed_get_phisc - - -subroutine ed_get_eimp(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(e_ineq))stop "ed_get_eimp error: e_ineq not allocated" - if(Nlat>size(e_ineq,1))stop "ed_get_eimp error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_eimp error: eimp has wrong rank' - rank (1) !all orbitals or all sites 1orbital - call assert_shape(self,[4],'ed_get_eimp','eimp') - self = [ed_Epot,ed_Eint,ed_Ehartree,ed_Eknot] - rank (2) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_eimp error: Nlat require for rank-2 array" - call assert_shape(self,[2,Nlat],'ed_get_eimp','eimp') - self = e_ineq - end select -end subroutine ed_get_eimp - - - - -subroutine ed_get_epot(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(e_ineq))stop "ed_get_epot error: e_ineq not allocated" - if(Nlat>size(e_ineq,1))stop "ed_get_epot error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_epot error: epot has wrong rank' - rank (0) - self = ed_Epot - rank (1) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_epot error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat],'ed_get_epot','epot') - self = e_ineq(:,1) - end select -end subroutine ed_get_epot - - -subroutine ed_get_eint(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(e_ineq))stop "ed_get_eint error: e_ineq not allocated" - if(Nlat>size(e_ineq,1))stop "ed_get_eint error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_eint error: eint has wrong rank' - rank (0) - self = ed_Eint - rank (1) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_eint error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat],'ed_get_eint','eint') - self = e_ineq(:,2) - end select -end subroutine ed_get_eint - -subroutine ed_get_ehartree(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(e_ineq))stop "ed_get_ehartree error: e_ineq not allocated" - if(Nlat>size(e_ineq,1))stop "ed_get_ehartree error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_ehartree error: ehartree has wrong rank' - rank (0) - self = ed_Ehartree - rank (1) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_ehartree error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat],'ed_get_ehartree','ehartree') - self = e_ineq(:,3) - end select -end subroutine ed_get_ehartree - - -subroutine ed_get_eknot(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(e_ineq))stop "ed_get_eknot error: e_ineq not allocated" - if(Nlat>size(e_ineq,1))stop "ed_get_eknot error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_eknot error: eknot has wrong rank' - rank (0) - self = ed_Eknot - rank (1) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_eknot error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat],'ed_get_eknot','eknot') - self = e_ineq(:,4) - end select -end subroutine ed_get_eknot - - - - - - - - - - -subroutine ed_get_doubles(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(dd_ineq))stop "ed_get_doubles error: dd_ineq not allocated" - if(Nlat>size(dd_ineq,1))stop "ed_get_doubles error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_doubles error: doubles has wrong rank' - rank (1) !all orbitals or all sites 1orbital - call assert_shape(self,[4],'ed_get_doubles','doubles') - self = [ed_Dust,ed_Dund,ed_Dse,ed_Dph] - rank (2) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_doubles error: Nlat require for rank-2 array" - call assert_shape(self,[2,Nlat],'ed_get_doubles','doubles') - self = dd_ineq - end select -end subroutine ed_get_doubles - -subroutine ed_get_dust(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(dd_ineq))stop "ed_get_dust error: dd_ineq not allocated" - if(Nlat>size(dd_ineq,1))stop "ed_get_dust error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_dust error: dust has wrong rank' - rank (0) - self = ed_Dust - rank (1) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_dust error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat],'ed_get_dust','dust') - self = dd_ineq(:,1) - end select -end subroutine ed_get_dust - - -subroutine ed_get_dund(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(dd_ineq))stop "ed_get_dund error: dd_ineq not allocated" - if(Nlat>size(dd_ineq,1))stop "ed_get_dund error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_dund error: dund has wrong rank' - rank (0) - self = ed_Dund - rank (1) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_dund error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat],'ed_get_dund','dund') - self = dd_ineq(:,2) - end select -end subroutine ed_get_dund - -subroutine ed_get_dse(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(dd_ineq))stop "ed_get_dse error: dd_ineq not allocated" - if(Nlat>size(dd_ineq,1))stop "ed_get_dse error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_dse error: dse has wrong rank' - rank (0) - self = ed_Dse - rank (1) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_dse error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat],'ed_get_dse','dse') - self = dd_ineq(:,3) - end select -end subroutine ed_get_dse - - -subroutine ed_get_dph(self,Nlat) - real(8),dimension(..) :: self - integer,optional :: Nlat - ! - if(present(Nlat))then - if(.not.allocated(dd_ineq))stop "ed_get_dph error: dd_ineq not allocated" - if(Nlat>size(dd_ineq,1))stop "ed_get_dph error: required N_sites > evaluated N_sites" - endif - select rank(self) - rank default;stop 'ed_get_dph error: dph has wrong rank' - rank (0) - self = ed_Dph - rank (1) !all orbitals and sites - if(.not.present(Nlat))stop "ed_get_dph error: Nlat require for rank-2 array" - call assert_shape(self,[Nlat],'ed_get_dph','dph') - self = dd_ineq(:,4) - end select -end subroutine ed_get_dph - - - - - -subroutine ed_get_neigen_total(nlii,Nlat) - integer :: Nlat - integer,dimension(Nlat) :: nlii - nlii=0d0 - if(allocated(neigen_total_ineq))then - if(Nlat>size(neigen_total_ineq)) stop "ed_get_docc error: required N_sites > evaluated N_sites" - nlii=neigen_total_ineq - endif -end subroutine ed_get_neigen_total - - diff --git a/src/ED_IO/get_phi.f90 b/src/ED_IO/get_phi.f90 new file mode 100644 index 00000000..c652a5e3 --- /dev/null +++ b/src/ED_IO/get_phi.f90 @@ -0,0 +1,37 @@ +subroutine ed_get_phisc_n0(self,iorb) + real(8) :: self + integer,optional :: iorb + integer :: iorb_ + iorb_=1;if(present(iorb))iorb_=iorb + if(iorb_>Norb)stop "ed_get_phisc error: orbital index > N_orbital" + self = ed_phisc(iorb_) +end subroutine ed_get_phisc_n0 + +subroutine ed_get_phisc_n1(self,iorb,Nlat) + real(8),dimension(:) :: self + integer,optional :: iorb,Nlat + integer :: iorb_ + iorb_=1;if(present(iorb))iorb_=iorb + if(iorb_>Norb)stop "ed_get_phisc error: orbital index > N_orbital" + if(present(Nlat))then + if(.not.allocated(phisc_ineq))stop "ed_get_phisc error: phisc_ineq not allocated" + if(Nlat>size(phisc_ineq,1))stop "ed_get_phisc error: required N_sites > evaluated N_sites" + endif + if(present(Nlat))then + call assert_shape(self,[Nlat],'ed_get_phisc','phisc') + self = phisc_ineq(:,iorb_) + else + call assert_shape(self,[Norb],'ed_get_phisc','phisc') + self = ed_phisc + endif +end subroutine ed_get_phisc_n1 + +subroutine ed_get_phisc_n2(self,Nlat) + real(8),dimension(:,:) :: self + integer :: Nlat + if(.not.allocated(phisc_ineq))stop "ed_get_phisc error: phisc_ineq not allocated" + if(Nlat>size(phisc_ineq,1))stop "ed_get_phisc error: required N_sites > evaluated N_sites" + call assert_shape(self,[Nlat,Norb],'ed_get_phisc','phisc') + self = phisc_ineq +end subroutine ed_get_phisc_n2 + diff --git a/src/ED_IO/get_sigma.f90 b/src/ED_IO/get_sigma.f90 index 27239f4a..4aed2b9a 100644 --- a/src/ED_IO/get_sigma.f90 +++ b/src/ED_IO/get_sigma.f90 @@ -1,9 +1,9 @@ -subroutine ed_get_sigma_site(self,axis,type) - complex(8),dimension(..),intent(inout) :: self - character(len=*),optional :: axis - character(len=*),optional :: type - character(len=1) :: axis_ - character(len=1) :: type_ +subroutine ed_get_sigma_site_n3(self,axis,type) + complex(8),dimension(:,:,:),intent(inout) :: self + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ axis_='m';if(present(axis))axis_=trim(axis) type_='n';if(present(type))type_=trim(type) select case(type_) @@ -12,57 +12,70 @@ subroutine ed_get_sigma_site(self,axis,type) select case(axis_) case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') - select rank(self) - rank default; stop "ed_get_sigma ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,Lmats],'ed_get_sigma','self') self = nn2so_reshape(impSmats,Nspin,Norb,Lmats) - rank (5) - call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lmats],'ed_get_sigma','self') - self = impSmats - end select case('r','R') - select rank(self) - rank default; stop "ed_get_sigma ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,Lreal],'ed_get_sigma','self') self = nn2so_reshape(impSreal,Nspin,Norb,Lmats) - rank (5) - call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lreal],'ed_get_sigma','self') - self = impSreal - end select end select case('a','A') select case(axis_) case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') - select rank(self) - rank default; stop "ed_get_sigma ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,Lmats],'ed_get_sigma','self') self = nn2so_reshape(impSAmats,Nspin,Norb,Lmats) - rank (5) - call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lmats],'ed_get_sigma','self') - self = impSAmats - end select case('r','R') - select rank(self) - rank default; stop "ed_get_sigma ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,Lreal],'ed_get_sigma','self') self = nn2so_reshape(impSAreal,Nspin,Norb,Lmats) - rank (5) + end select + end select +end subroutine ed_get_sigma_site_n3 + + + +subroutine ed_get_sigma_site_n5(self,axis,type) + complex(8),dimension(:,:,:,:,:),intent(inout) :: self + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + select case(type_) + case default; stop "ed_get_sigma ERROR: type is neither Normal, nor Anomalous" + case ('n','N') + select case(axis_) + case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lmats],'ed_get_sigma','self') + self = impSmats + case('r','R') + call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lreal],'ed_get_sigma','self') + self = impSreal + end select + case('a','A') + select case(axis_) + case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lmats],'ed_get_sigma','self') + self = impSAmats + case('r','R') call assert_shape(self,[Nspin,Nspin,Norb,Norb,Lreal],'ed_get_sigma','self') self = impSAreal - end select end select end select -end subroutine ed_get_sigma_site +end subroutine ed_get_sigma_site_n5 + + +!################################################################## +!################################################################## +!################################################################## -subroutine ed_get_sigma_lattice(self,nlat,axis,type) - complex(8),dimension(..),intent(inout) :: self + +subroutine ed_get_sigma_lattice_n3(self,nlat,axis,type) + complex(8),dimension(:,:,:),intent(inout) :: self integer,intent(in) :: nlat character(len=*),optional :: axis character(len=*),optional :: type @@ -78,72 +91,108 @@ subroutine ed_get_sigma_lattice(self,nlat,axis,type) select case(axis_) case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') - select rank(self) - rank default; stop "ed_get_sigma ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lmats],'ed_get_sigma','self') self = nnn2lso_reshape(Smats_ineq,Nlat,Nspin,Norb,Lmats) - rank (4) + case('r','R') + call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lreal],'ed_get_sigma','self') + self = nnn2lso_reshape(Sreal_ineq,Nlat,Nspin,Norb,Lreal) + end select + case('a','A') + select case(axis_) + case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lmats],'ed_get_sigma','self') + self = nnn2lso_reshape(SAmats_ineq,Nlat,Nspin,Norb,Lmats) + case('r','R') + call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lreal],'ed_get_sigma','self') + self = nnn2lso_reshape(SAreal_ineq,Nlat,Nspin,Norb,Lreal) + end select + end select +end subroutine ed_get_sigma_lattice_n3 + + +subroutine ed_get_sigma_lattice_n4(self,nlat,axis,type) + complex(8),dimension(:,:,:,:),intent(inout) :: self + integer,intent(in) :: nlat + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + integer :: ilat + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + if(Nlat/=size(Smats_ineq,1))stop "ERROR ed_get_sigma: wrong Nlat" + select case(type_) + case default; stop "ed_get_sigma ERROR: type is neither Normal, nor Anomalous" + case ('n','N') + select case(axis_) + case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') call assert_shape(self,[Nlat,Nspin*Norb,Nspin*Norb,Lmats],'ed_get_sigma','self') do ilat=1,Nlat self(ilat,:,:,:) = nn2so_reshape(Smats_ineq(ilat,:,:,:,:,:),Nspin,Norb,Lmats) enddo - rank (6) - call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lmats],'ed_get_sigma','self') - self = Smats_ineq - end select case('r','R') - select rank(self) - rank default; stop "ed_get_sigma ERROR: self has a wrong rank" - rank (3) - call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lreal],'ed_get_sigma','self') - self = nnn2lso_reshape(Sreal_ineq,Nlat,Nspin,Norb,Lreal) - rank (4) call assert_shape(self,[Nlat,Nspin*Norb,Nspin*Norb,Lreal],'ed_get_sigma','self') do ilat=1,Nlat self(ilat,:,:,:) = nn2so_reshape(Sreal_ineq(ilat,:,:,:,:,:),Nspin,Norb,Lreal) enddo - rank (6) - call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lreal],'ed_get_sigma','self') - self = Sreal_ineq - end select end select case('a','A') select case(axis_) case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" case ('m','M') - select rank(self) - rank default; stop "ed_get_sigma ERROR: self has a wrong rank" - rank (3) - call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lmats],'ed_get_sigma','self') - self = nnn2lso_reshape(SAmats_ineq,Nlat,Nspin,Norb,Lmats) - rank (4) call assert_shape(self,[Nlat,Nspin*Norb,Nspin*Norb,Lmats],'ed_get_sigma','self') do ilat=1,Nlat self(ilat,:,:,:) = nn2so_reshape(SAmats_ineq(ilat,:,:,:,:,:),Nspin,Norb,Lmats) enddo - rank (6) - call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lmats],'ed_get_sigma','self') - self = SAmats_ineq - end select case('r','R') - select rank(self) - rank default; stop "ed_get_sigma ERROR: self has a wrong rank" - rank (3) - call assert_shape(self,[Nlat*Nspin*Norb,Nlat*Nspin*Norb,Lreal],'ed_get_sigma','self') - self = nnn2lso_reshape(SAreal_ineq,Nlat,Nspin,Norb,Lreal) - rank (4) call assert_shape(self,[Nlat,Nspin*Norb,Nspin*Norb,Lreal],'ed_get_sigma','self') do ilat=1,Nlat self(ilat,:,:,:) = nn2so_reshape(SAreal_ineq(ilat,:,:,:,:,:),Nspin,Norb,Lreal) enddo - rank (6) + end select + end select +end subroutine ed_get_sigma_lattice_n4 + + + + + +subroutine ed_get_sigma_lattice_n6(self,nlat,axis,type) + complex(8),dimension(:,:,:,:,:,:),intent(inout) :: self + integer,intent(in) :: nlat + character(len=*),optional :: axis + character(len=*),optional :: type + character(len=1) :: axis_ + character(len=1) :: type_ + integer :: ilat + axis_='m';if(present(axis))axis_=trim(axis) + type_='n';if(present(type))type_=trim(type) + if(Nlat/=size(Smats_ineq,1))stop "ERROR ed_get_sigma: wrong Nlat" + select case(type_) + case default; stop "ed_get_sigma ERROR: type is neither Normal, nor Anomalous" + case ('n','N') + select case(axis_) + case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lmats],'ed_get_sigma','self') + self = Smats_ineq + case('r','R') + call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lreal],'ed_get_sigma','self') + self = Sreal_ineq + end select + case('a','A') + select case(axis_) + case default;stop "ed_get_sigma ERROR: axis is neither Matsubara, nor Realaxis" + case ('m','M') + call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lmats],'ed_get_sigma','self') + self = SAmats_ineq + case('r','R') call assert_shape(self,[Nlat,Nspin,Nspin,Norb,Norb,Lreal],'ed_get_sigma','self') self = SAreal_ineq - end select end select end select -end subroutine ed_get_sigma_lattice - +end subroutine ed_get_sigma_lattice_n6 diff --git a/src/ED_IO/rebuild_impG.f90 b/src/ED_IO/rebuild_impG.f90 index ea74f940..c45beee0 100644 --- a/src/ED_IO/rebuild_impG.f90 +++ b/src/ED_IO/rebuild_impG.f90 @@ -1,7 +1,7 @@ -subroutine rebuild_gimp_single(zeta,gimp,fimp) +subroutine rebuild_gimp_single_n3(zeta,gimp,fimp) complex(8),dimension(:) :: zeta - complex(8),dimension(..) :: Gimp - complex(8),dimension(..),optional :: Fimp + complex(8),dimension(:,:,:) :: Gimp + complex(8),dimension(:,:,:),optional :: Fimp integer :: i,L logical :: check complex(8),dimension(:,:,:,:,:),allocatable :: gf,ff @@ -24,36 +24,61 @@ subroutine rebuild_gimp_single(zeta,gimp,fimp) ! call deallocate_GFmatrix(impGmatrix) ! - select rank(gimp) - rank default; stop "rebuild_gimp ERROR: gimp has a wrong rank" - rank (3) call assert_shape(gimp,[Nspin*Norb,Nspin*Norb,L],'rebuild_gimp','gimp') gimp = nn2so_reshape(gf,Nspin,Norb,L) - rank (5) - call assert_shape(gimp,[Nspin,Nspin,Norb,Norb,L],'rebuild_gimp','gimp') - gimp = gf - end select if(present(fimp))then - select rank(fimp) - rank default; stop "rebuild_gimp ERROR: fimp has a wrong rank" - rank (3) call assert_shape(fimp,[Nspin*Norb,Nspin*Norb,L],'rebuild_gimp','fimp') fimp = nn2so_reshape(ff,Nspin,Norb,L) - rank (5) + endif + return +end subroutine rebuild_gimp_single_n3 + +subroutine rebuild_gimp_single_n5(zeta,gimp,fimp) + complex(8),dimension(:) :: zeta + complex(8),dimension(:,:,:,:,:) :: Gimp + complex(8),dimension(:,:,:,:,:),optional :: Fimp + integer :: i,L + logical :: check + complex(8),dimension(:,:,:,:,:),allocatable :: gf,ff + ! + call ed_read_impGmatrix() + ! + L = size(zeta) + allocate(gf(Nspin,Nspin,Norb,Norb,L)) + allocate(ff(Nspin,Nspin,Norb,Norb,L)) + ! + select case(ed_mode) + case default; + call rebuild_gimp_normal(zeta,gf) + case("superc"); + if(.not.present(fimp))stop "rebuild_gimp ERROR: fimp not present in ed_mode=superc" + call rebuild_gimp_superc(zeta,gf,ff) + case("nonsu2"); + call rebuild_gimp_nonsu2(zeta,gf) + end select + ! + call deallocate_GFmatrix(impGmatrix) + ! + call assert_shape(gimp,[Nspin,Nspin,Norb,Norb,L],'rebuild_gimp','gimp') + gimp = gf + if(present(fimp))then call assert_shape(fimp,[Nspin,Nspin,Norb,Norb,L],'rebuild_gimp','fimp') fimp = ff - end select endif return -end subroutine rebuild_gimp_single +end subroutine rebuild_gimp_single_n5 +!################################################################## +!################################################################## +!################################################################## -subroutine rebuild_gimp_ineq(zeta,Nlat,gimp,fimp) + +subroutine rebuild_gimp_ineq_n3(zeta,Nlat,gimp,fimp) complex(8),dimension(:) :: zeta integer :: Nlat - complex(8),dimension(..) :: Gimp - complex(8),dimension(..),optional :: Fimp + complex(8),dimension(:,:,:) :: Gimp + complex(8),dimension(:,:,:),optional :: Fimp integer :: i,ilat,Nineq,L logical :: check complex(8),dimension(:,:,:,:,:,:),allocatable :: gf,ff @@ -67,49 +92,104 @@ subroutine rebuild_gimp_ineq(zeta,Nlat,gimp,fimp) do ilat = 1, Nineq ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) if(present(fimp))then - call rebuild_gimp_single(zeta,& + call rebuild_gimp_single_n5(zeta,& gf(ilat,:,:,:,:,:),& ff(ilat,:,:,:,:,:)) else - call rebuild_gimp_single(zeta,& + call rebuild_gimp_single_n5(zeta,& ff(ilat,:,:,:,:,:)) endif enddo call ed_reset_suffix ! - select rank(gimp) - rank default; stop "rebuild_gimp ERROR: gimp has a wrong rank" - rank (3) call assert_shape(gimp,[Nineq*Nspin*Norb,Nineq*Nspin*Norb,L],'rebuild_gimp','gimp') gimp = nnn2lso_reshape(gf,Nineq,Nspin,Norb,L) - rank (4) + + if(present(fimp))then + call assert_shape(fimp,[Nineq*Nspin*Norb,Nineq*Nspin*Norb,L],'rebuild_gimp','fimp') + fimp = nnn2lso_reshape(ff,Nineq,Nspin,Norb,L) + endif + return +end subroutine rebuild_gimp_ineq_n3 + +subroutine rebuild_gimp_ineq_n4(zeta,Nlat,gimp,fimp) + complex(8),dimension(:) :: zeta + integer :: Nlat + complex(8),dimension(:,:,:,:) :: Gimp + complex(8),dimension(:,:,:,:),optional :: Fimp + integer :: i,ilat,Nineq,L + logical :: check + complex(8),dimension(:,:,:,:,:,:),allocatable :: gf,ff + ! + Nineq=Nlat + L=size(zeta) + ! + allocate(Gf(Nineq,Nspin,Nspin,Norb,Norb,L)) + allocate(Ff(Nineq,Nspin,Nspin,Norb,Norb,L)) + ! + do ilat = 1, Nineq + ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) + if(present(fimp))then + call rebuild_gimp_single_n5(zeta,& + gf(ilat,:,:,:,:,:),& + ff(ilat,:,:,:,:,:)) + else + call rebuild_gimp_single_n5(zeta,& + ff(ilat,:,:,:,:,:)) + endif + enddo + call ed_reset_suffix + ! call assert_shape(gimp,[Nineq,Nspin*Norb,Nspin*Norb,L],'rebuild_gimp','gimp') do ilat=1,Nineq gimp(ilat,:,:,:) = nn2so_reshape(gf(ilat,:,:,:,:,:),Nspin,Norb,L) enddo - rank (6) - call assert_shape(gimp,[Nineq,Nspin,Nspin,Norb,Norb,L],'rebuild_gimp','gimp') - gimp = gf - end select if(present(fimp))then - select rank(fimp) - rank default; stop "rebuild_gimp ERROR: fimp has a wrong rank" - rank (3) - call assert_shape(fimp,[Nineq*Nspin*Norb,Nineq*Nspin*Norb,L],'rebuild_gimp','fimp') - fimp = nnn2lso_reshape(ff,Nineq,Nspin,Norb,L) - rank (4) call assert_shape(fimp,[Nineq,Nspin*Norb,Nspin*Norb,L],'rebuild_gimp','fimp') do ilat=1,Nineq fimp(ilat,:,:,:) = nn2so_reshape(ff(ilat,:,:,:,:,:),Nspin,Norb,L) enddo - rank (6) + endif + return +end subroutine rebuild_gimp_ineq_n4 + + +subroutine rebuild_gimp_ineq_n6(zeta,Nlat,gimp,fimp) + complex(8),dimension(:) :: zeta + integer :: Nlat + complex(8),dimension(:,:,:,:,:,:) :: Gimp + complex(8),dimension(:,:,:,:,:,:),optional :: Fimp + integer :: i,ilat,Nineq,L + logical :: check + complex(8),dimension(:,:,:,:,:,:),allocatable :: gf,ff + ! + Nineq=Nlat + L=size(zeta) + ! + allocate(Gf(Nineq,Nspin,Nspin,Norb,Norb,L)) + allocate(Ff(Nineq,Nspin,Nspin,Norb,Norb,L)) + ! + do ilat = 1, Nineq + ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) + if(present(fimp))then + call rebuild_gimp_single_n5(zeta,& + gf(ilat,:,:,:,:,:),& + ff(ilat,:,:,:,:,:)) + else + call rebuild_gimp_single_n5(zeta,& + ff(ilat,:,:,:,:,:)) + endif + enddo + call ed_reset_suffix + ! + call assert_shape(gimp,[Nineq,Nspin,Nspin,Norb,Norb,L],'rebuild_gimp','gimp') + gimp = gf + if(present(fimp))then call assert_shape(fimp,[Nineq,Nspin,Nspin,Norb,Norb,L],'rebuild_gimp','fimp') fimp = ff - end select endif return -end subroutine rebuild_gimp_ineq - +end subroutine rebuild_gimp_ineq_n6 diff --git a/src/ED_IO/rebuild_impSigma.f90 b/src/ED_IO/rebuild_impSigma.f90 index bd7b8624..2ba37cb0 100644 --- a/src/ED_IO/rebuild_impSigma.f90 +++ b/src/ED_IO/rebuild_impSigma.f90 @@ -1,7 +1,7 @@ -subroutine rebuild_sigma_single(zeta,sigma,self) +subroutine rebuild_sigma_single_n3(zeta,sigma,self) complex(8),dimension(:) :: zeta - complex(8),dimension(..) :: Sigma - complex(8),dimension(..),optional :: Self + complex(8),dimension(:,:,:) :: Sigma + complex(8),dimension(:,:,:),optional :: Self integer :: i,L logical :: check complex(8),dimension(:,:,:,:,:),allocatable :: sn,sa @@ -29,36 +29,67 @@ subroutine rebuild_sigma_single(zeta,sigma,self) call deallocate_dmft_bath(dmft_bath) call deallocate_GFmatrix(impGmatrix) ! - select rank(sigma) - rank default; stop "rebuild_sigma ERROR: sigma has a wrong rank" - rank (3) call assert_shape(sigma,[Nspin*Norb,Nspin*Norb,L],'rebuild_sigma','sigma') sigma = nn2so_reshape(sn,Nspin,Norb,L) - rank (5) - call assert_shape(sigma,[Nspin,Nspin,Norb,Norb,L],'rebuild_sigma','sigma') - sigma = sn - end select if(present(self))then - select rank(self) - rank default; stop "rebuild_sigma ERROR: self has a wrong rank" - rank (3) call assert_shape(self,[Nspin*Norb,Nspin*Norb,L],'rebuild_sigma','self') self = nn2so_reshape(sa,Nspin,Norb,L) - rank (5) + endif + return +end subroutine rebuild_sigma_single_n3 + +subroutine rebuild_sigma_single_n5(zeta,sigma,self) + complex(8),dimension(:) :: zeta + complex(8),dimension(:,:,:,:,:) :: Sigma + complex(8),dimension(:,:,:,:,:),optional :: Self + integer :: i,L + logical :: check + complex(8),dimension(:,:,:,:,:),allocatable :: sn,sa + ! + call ed_read_impGmatrix() + ! + if(.not.allocated(impHloc))stop "rebuild_sigma error: impHloc not allocated. Call ed_set_Hloc first." + ! + call allocate_dmft_bath(dmft_bath) + call init_dmft_bath(dmft_bath,used=.true.) + ! + L = size(zeta) + allocate(sn(Nspin,Nspin,Norb,Norb,L)) + allocate(sa(Nspin,Nspin,Norb,Norb,L)) + select case(ed_mode) + case default; + call rebuild_sigma_normal(zeta,sn) + case("superc"); + if(.not.present(self))stop "rebuild_sigma ERROR: self not present in ed_mode=superc" + call rebuild_sigma_superc(zeta,sn,sa) + case("nonsu2"); + call rebuild_sigma_nonsu2(zeta,sn) + end select + ! + call deallocate_dmft_bath(dmft_bath) + call deallocate_GFmatrix(impGmatrix) + ! + call assert_shape(sigma,[Nspin,Nspin,Norb,Norb,L],'rebuild_sigma','sigma') + sigma = sn + if(present(self))then call assert_shape(self,[Nspin,Nspin,Norb,Norb,L],'rebuild_sigma','self') self = sa - end select endif return -end subroutine rebuild_sigma_single +end subroutine rebuild_sigma_single_n5 + + +!################################################################## +!################################################################## +!################################################################## -subroutine rebuild_sigma_ineq(zeta,Nlat,sigma,self) +subroutine rebuild_sigma_ineq_n3(zeta,Nlat,sigma,self) complex(8),dimension(:) :: zeta integer :: Nlat - complex(8),dimension(..) :: Sigma - complex(8),dimension(..),optional :: Self + complex(8),dimension(:,:,:) :: Sigma + complex(8),dimension(:,:,:),optional :: Self integer :: i,ilat,Nineq,L logical :: check complex(8),dimension(:,:,:,:,:,:),allocatable :: sn,sa @@ -75,46 +106,100 @@ subroutine rebuild_sigma_ineq(zeta,Nlat,sigma,self) ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) call ed_set_Hloc(Hloc_ineq(ilat,:,:,:,:)) if(present(self))then - call rebuild_sigma_single(zeta,sn(ilat,:,:,:,:,:),sa(ilat,:,:,:,:,:)) + call rebuild_sigma_single_n5(zeta,sn(ilat,:,:,:,:,:),sa(ilat,:,:,:,:,:)) else - call rebuild_sigma_single(zeta,sa(ilat,:,:,:,:,:)) + call rebuild_sigma_single_n5(zeta,sa(ilat,:,:,:,:,:)) endif enddo call ed_reset_suffix ! - select rank(sigma) - rank default; stop "rebuild_sigma ERROR: sigma has a wrong rank" - rank (3) call assert_shape(sigma,[Nineq*Nspin*Norb,Nineq*Nspin*Norb,L],'rebuild_sigma','sigma') sigma = nnn2lso_reshape(sn,Nineq,Nspin,Norb,L) - rank (4) + if(present(self))then + call assert_shape(self,[Nineq*Nspin*Norb,Nineq*Nspin*Norb,L],'rebuild_sigma','self') + self = nnn2lso_reshape(sa,Nineq,Nspin,Norb,L) + endif + return +end subroutine rebuild_sigma_ineq_n3 + +subroutine rebuild_sigma_ineq_n4(zeta,Nlat,sigma,self) + complex(8),dimension(:) :: zeta + integer :: Nlat + complex(8),dimension(:,:,:,:) :: Sigma + complex(8),dimension(:,:,:,:),optional :: Self + integer :: i,ilat,Nineq,L + logical :: check + complex(8),dimension(:,:,:,:,:,:),allocatable :: sn,sa + ! + Nineq=Nlat + L=size(zeta) + ! + if(.not.allocated(Hloc_ineq))stop "ed_rebuild_gf error: Hloc_ineq not allocated. Call ed_set_Hloc first." + ! + allocate(Sn(Nineq,Nspin,Nspin,Norb,Norb,L)) + allocate(Sa(Nineq,Nspin,Nspin,Norb,Norb,L)) + ! + do ilat = 1, Nineq + ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) + call ed_set_Hloc(Hloc_ineq(ilat,:,:,:,:)) + if(present(self))then + call rebuild_sigma_single_n5(zeta,sn(ilat,:,:,:,:,:),sa(ilat,:,:,:,:,:)) + else + call rebuild_sigma_single_n5(zeta,sa(ilat,:,:,:,:,:)) + endif + enddo + call ed_reset_suffix + ! call assert_shape(sigma,[Nineq,Nspin*Norb,Nspin*Norb,L],'rebuild_sigma','sigma') do ilat=1,Nineq sigma(ilat,:,:,:) = nn2so_reshape(sn(ilat,:,:,:,:,:),Nspin,Norb,L) enddo - rank (6) - call assert_shape(sigma,[Nineq,Nspin,Nspin,Norb,Norb,L],'rebuild_sigma','sigma') - sigma = sn - end select if(present(self))then - select rank(self) - rank default; stop "rebuild_sigma ERROR: self has a wrong rank" - rank (3) - call assert_shape(self,[Nineq*Nspin*Norb,Nineq*Nspin*Norb,L],'rebuild_sigma','self') - self = nnn2lso_reshape(sa,Nineq,Nspin,Norb,L) - rank (4) call assert_shape(self,[Nineq,Nspin*Norb,Nspin*Norb,L],'rebuild_sigma','self') do ilat=1,Nineq self(ilat,:,:,:) = nn2so_reshape(sa(ilat,:,:,:,:,:),Nspin,Norb,L) enddo - rank (6) + endif + return +end subroutine rebuild_sigma_ineq_n4 + + +subroutine rebuild_sigma_ineq_n6(zeta,Nlat,sigma,self) + complex(8),dimension(:) :: zeta + integer :: Nlat + complex(8),dimension(:,:,:,:,:,:) :: Sigma + complex(8),dimension(:,:,:,:,:,:),optional :: Self + integer :: i,ilat,Nineq,L + logical :: check + complex(8),dimension(:,:,:,:,:,:),allocatable :: sn,sa + ! + Nineq=Nlat + L=size(zeta) + ! + if(.not.allocated(Hloc_ineq))stop "ed_rebuild_gf error: Hloc_ineq not allocated. Call ed_set_Hloc first." + ! + allocate(Sn(Nineq,Nspin,Nspin,Norb,Norb,L)) + allocate(Sa(Nineq,Nspin,Nspin,Norb,Norb,L)) + ! + do ilat = 1, Nineq + ed_file_suffix=reg(ineq_site_suffix)//str(ilat,site_indx_padding) + call ed_set_Hloc(Hloc_ineq(ilat,:,:,:,:)) + if(present(self))then + call rebuild_sigma_single_n5(zeta,sn(ilat,:,:,:,:,:),sa(ilat,:,:,:,:,:)) + else + call rebuild_sigma_single_n5(zeta,sa(ilat,:,:,:,:,:)) + endif + enddo + call ed_reset_suffix + ! + call assert_shape(sigma,[Nineq,Nspin,Nspin,Norb,Norb,L],'rebuild_sigma','sigma') + sigma = sn + if(present(self))then call assert_shape(self,[Nineq,Nspin,Nspin,Norb,Norb,L],'rebuild_sigma','self') self = sa - end select endif return -end subroutine rebuild_sigma_ineq - +end subroutine rebuild_sigma_ineq_n6