Skip to content

Commit

Permalink
Intermediate commit
Browse files Browse the repository at this point in the history 
- code cleanup
- update tests: by using numpy broadcast when calculating Gloc, python
  is almost as fast as fortran
- update tests: added Fortran example. The makefile doesn't work
  • Loading branch information
@lcrippa
lcrippa committed Nov 8, 2024
1 parent 55bbbe1 commit f4be00f
Show file tree
Hide file tree
Showing 6 changed files with 297 additions and 24 deletions.
14 changes: 2 additions & 12 deletions src/python/edipy2/func_aux_funx.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -207,28 +207,18 @@ def check_convergence(self,func,threshold,N1=None,N2=None):
colorprefix= self.BOLD + self.YELLOW
else:
colorprefix= self.BOLD + self.RED

print("\n")

if self.whichiter <= N2:
if np.prod(np.shape(errvec)) > 1:
print(colorprefix + "max error=" + self.COLOREND + f"{errmax:.6e}")
print(colorprefix + " "*(np.prod(np.shape(errvec)) > 1)+"error=" + self.COLOREND + f"{err:.6e}")
if np.prod(np.shape(errvec)) > 1:
print(colorprefix + "min error=" + self.COLOREND + f"{errmin:.6e}\n")
print(colorprefix + "min error=" + self.COLOREND + f"{errmin:.6e}")
else:
print("Not converged after "+str(N2)+" iterations.")
with open("ERROR.ERR", "w") as file:
file.write("Not converged after "+str(N2)+" iterations.")
print("\n")

#delete methods after convergence or Nmax exceeded
#Disabled: I think these need to remain defined in case of fixed-density calculations where
#the truth value of the convergence flag is modified afterwards.
#Moving this to finalize_environment
#if conv_bool:
# del self.oldfunc
# del self.gooditer
# del self.whichiter

#pass to other cores:
try:
Expand Down
17 changes: 17 additions & 0 deletions test/fortran/Makefile.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
FC=mpif90

INC=$(shell pkg-config --cflags scifor edipack2)
LIB=$(shell pkg-config --libs scifor edipack2)

FFLAG = -O2 -ffree-line-length-none -cpp -D_MPI

# Extends the implicit support of the Makefile to .f90 files
.SUFFIXES: .f90

all:
@echo ""
$(FC) $(FFLAG) $(INC) $(LIB) hm_bethe.f90 -o hm_bethe.exe

clean:
@echo "Cleaning:"
@rm -f *.exe *.mod *.o *~
8 changes: 8 additions & 0 deletions test/fortran/hamiltonian.restart
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
#Ek_l1_s1 Vk_l1_s1
-1.081558684145E+00 2.909966891416E-01
-1.309724255302E-01 -1.697488530953E-01
-1.815378548632E-02 8.126987899873E-02
1.838474629692E-09 3.390315786474E-02
1.815378762971E-02 8.126985492899E-02
1.309723274470E-01 -1.697487952453E-01
1.081558659122E+00 2.909966975902E-01
189 changes: 189 additions & 0 deletions test/fortran/hm_bethe.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,189 @@
program hm_bethe
USE EDIPACK2
USE SCIFOR
USE MPI
implicit none
!INPUT:
integer :: Le
real(8) :: wmixing
real(8) :: Wband
!
!Bath:
integer :: Nb
real(8),dimension(:),allocatable :: Bath,Bath_prev
!Local Variables:
complex(8),dimension(:,:,:,:),allocatable :: Hloc
real(8),dimension(:),allocatable :: Ebands
real(8),dimension(:),allocatable :: Dbands

real(8) :: de
!Local Functions:
complex(8),allocatable,dimension(:,:,:,:,:) :: Smats,Sreal
complex(8),allocatable,dimension(:,:,:,:,:) :: Gmats,Greal
complex(8),allocatable,dimension(:,:,:,:,:) :: Delta
!Other:
integer :: iloop,iorb,i
logical :: converged
complex(8) :: zeta
real(8),dimension(:),allocatable :: wm,wr
integer :: comm,rank
logical :: master


!MPI Wrapper from SciFortran:
call init_MPI()
comm = MPI_COMM_WORLD
call StartMsg_MPI(comm)
rank = get_Rank_MPI(comm)
master = get_Master_MPI(comm)


!READ ED INPUT:
call ed_read_input('inputED.conf')
if(Nspin/=1.OR.Norb/=1)stop "This test code is for Nspin=1 + Norb=1."
Le = 1000
wmixing= 0.5d0
Wband = 1d0

!BUILD Density of States:
!linspace, dens_bethe, arange, pi in SciFortran
allocate(Ebands(Le))
allocate(Dbands(Le))
Ebands = linspace(-Wband,Wband,Le,mesh=de)
Dbands = dens_bethe(Ebands,Wband)*de

!BUILD frequency array:
allocate(wm(Lmats),wr(Lreal))
wm = pi/beta*(2*arange(1,Lmats)-1)
wr = linspace(wini,wfin,Lreal)


!ALL functions must have shape [Nspin,Nspin,Norb,Norb(,L)]:
allocate(Delta(Nspin,Nspin,Norb,Norb,Lmats))
allocate(Smats(Nspin,Nspin,Norb,Norb,Lmats))
allocate(Gmats(Nspin,Nspin,Norb,Norb,Lmats))
allocate(Sreal(Nspin,Nspin,Norb,Norb,Lreal))
allocate(Greal(Nspin,Nspin,Norb,Norb,Lreal))
allocate(Hloc(Nspin,Nspin,Norb,Norb))
Hloc = 0d0


!SETUP SOLVER
Nb=ed_get_bath_dimension()
allocate(bath(Nb),bath_prev(Nb))
call ed_set_hloc(Hloc)
call ed_init_solver(bath)

!DMFT CYCLE
iloop=0;converged=.false.
do while(.not.converged.AND.iloop<nloop)
iloop=iloop+1
if(master)write(*,*)"DMFT-loop:",iloop,"/",nloop
!
!Solve the EFFECTIVE IMPURITY PROBLEM (first w/ a guess for the bath)
call ed_solve(bath)
!
!Get Self-energies
call ed_get_sigma(Smats,axis="m")
call ed_get_sigma(Sreal,axis="r")
!
!Compute the local gf:
do i=1,Lmats
zeta = xi*wm(i) + xmu - Smats(1,1,1,1,i)
Gmats(1,1,1,1,i) = sum(Dbands(:)/(zeta - Ebands(:)))
enddo
do i=1,Lreal
zeta = cmplx(wr(i),eps) + xmu - Sreal(1,1,1,1,i)
Greal(1,1,1,1,i) = sum(Dbands(:)/(zeta - Ebands(:)))
enddo
if(master)call splot("Gloc_iw.dat",wm,Gmats(1,1,1,1,:)) !splot is in SciFortran
if(master)call splot("Gloc_realw.dat",wr,Greal(1,1,1,1,:)) !splot is in SciFortran
!
!Get the Delta function and FIT:
Delta(1,1,1,1,:) = 0.25d0*Wband*Gmats(1,1,1,1,:)
call ed_chi2_fitgf(Delta,Bath,ispin=1)
!
!MIXING:
if(iloop>1)Bath = wmixing*Bath + (1.d0-wmixing)*Bath_prev
Bath_prev=Bath
!
!Check convergence
converged = check_convergence(Delta(1,1,1,1,:),dmft_error,nsuccess,nloop)
if(master)write(*,*)"---------"
enddo

call finalize_MPI()


contains

function check_convergence(Xnew,eps,N1,N2) result(convergence)
complex(8),intent(in) :: Xnew(:)
real(8),intent(in) :: eps
integer,intent(in) :: N1,N2
logical :: convergence
integer :: i,j,Msum
real(8) :: err
real(8) :: M,S
complex(8),save,allocatable :: Xold(:)
integer,save :: success=0,check=1
character(len=100) :: file_
file_='error.err'
Msum=size(Xnew)
if(.not.allocated(Xold))then
allocate(Xold(Msum))
Xold=0.d0
endif
S=0.d0 ; M=0.d0
do i=1,Msum
M=M + abs(Xnew(i)-Xold(i))
S=S + abs(Xnew(i))
enddo
err= M/S
Xold=Xnew
if(master)then
open(10,file=reg(file_),position="append")
write(10,"(I5,ES15.7)")check,err
close(10)
endif
if(err < eps)then
success=success+1
else
success=0
endif
convergence=.false.
if(success > N1)convergence=.true.
if(check>=N2)then
if(master)then
open(10,file="ERROR.README")
write(10,*)""
close(10)
write(*,"(A,I4,A)")"Not converged after",N2," iterations."
endif
endif
if(convergence)then
if(master)write(*,"(A,ES15.7)")bold_green("error="),err
else
if(err < eps)then
if(master)write(*,"(A,ES15.7)")bold_yellow("error="),err
else
if(master)write(*,"(A,ES15.7)")bold_red("error="),err
endif
endif
check=check+1
end function check_convergence


end program hm_bethe












74 changes: 74 additions & 0 deletions test/fortran/inputED.conf
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
NORB=1 !Number of impurity orbitals (max 5).
NBATH=7 !Number of bath sites:(normal=>Nbath per orb)(hybrid=>Nbath total)(replica=>Nbath=Nreplica)
NSPIN=1 !Number of spin degeneracy (max 2)
PH_TYPE=1 !Shape e-ph interaction: 1=orbital density, 2=orbital hybridization
NPH=0 !Max number of phonons allowed (cut off)
ULOC=2.000000000,0.d0,0.d0,0.d0,0.d0 !Values of the local interaction per orbital (max 5)
UST=0d0 !Value of the inter-orbital interaction term
JH=0.d0 !Hunds coupling
JX=0.d0 !S-E coupling
JP=0.d0 !P-H coupling
BETA=1000.000000000 !Inverse temperature, at T=0 is used as a IR cut-off.
XMU=0.d0 !Chemical potential. If HFMODE=T, xmu=0 indicates half-filling condition.
NLOOP=1 !Max number of DMFT iterations.
DMFT_ERROR=1.000000000E-04 !Error threshold for DMFT convergence
SB_FIELD=1.000000000E-01 !Value of a symmetry breaking field for magnetic solutions.
ED_DIAG_TYPE=lanc !flag to select the diagonalization type: 'lanc' for Lanczos/Davidson, 'full' for Full diagonalization method
ED_FINITE_TEMP=F !flag to select finite temperature method. note that if T then lanc_nstates_total must be > 1
ED_TWIN=T !flag to reduce (T) or not (F,default) the number of visited sector using twin symmetry.
ED_SECTORS=F !flag to reduce sector scan for the spectrum to specific sectors +/- ed_sectors_shift.
ED_SECTORS_SHIFT=1 !shift to ed_sectors
ED_SPARSE_H=T !flag to select storage of sparse matrix H (mem--, cpu++) if TRUE, or direct on-the-fly H*v product (mem++, cpu--) if FALSE
ED_TOTAL_UD=T !flag to select which type of quantum numbers have to be considered: T (default) total Nup-Ndw, F orbital based Nup-Ndw
ED_SOLVE_OFFDIAG_GF=F !flag to select the calculation of the off-diagonal impurity GF. this is T by default if bath_type/=normal
ED_PRINT_SIGMA=T !flag to print impurity Self-energies
ED_PRINT_G=T !flag to print impurity Greens function
ED_PRINT_G0=F !flag to print non-interacting impurity Greens function
ED_VERBOSE=4 !Verbosity level: 0=almost nothing --> 5:all. Really: all
NSUCCESS=1 !Number of successive iterations below threshold for convergence
LMATS=6000 !Number of Matsubara frequencies.
LREAL=4000 !Number of real-axis frequencies.
LTAU=1000 !Number of imaginary time points.
LFIT=1000 !Number of Matsubara frequencies used in the \Chi2 fit.
LPOS=100 !Number of points for the lattice PDF.
NREAD=0.d0 !Objective density for fixed density calculations.
NERR=1.000000000E-04 !Error threshold for fixed density calculations.
NDELTA=1.000000000E-01 !Initial step for fixed density calculations.
NCOEFF=1.000000000 !multiplier for the initial ndelta read from a file (ndelta-->ndelta*ncoeff).
WINI=-5.000000000 !Smallest real-axis frequency
WFIN=5.000000000 !Largest real-axis frequency
XMIN=-3.000000000 !Smallest position for the lattice PDF
XMAX=3.000000000 !Largest position for the lattice PDF
CHISPIN_FLAG=F !Flag to activate spin susceptibility calculation.
CHIDENS_FLAG=F !Flag to activate density susceptibility calculation.
CHIPAIR_FLAG=F !Flag to activate pair susceptibility calculation.
CHIEXCT_FLAG=F !Flag to activate excitonis susceptibility calculation.
HFMODE=T !Flag to set the Hartree form of the interaction (n-1/2). see xmu.
EPS=4.000000000E-02 !Broadening on the real-axis.
CUTOFF=1.000000000E-09 !Spectrum cut-off, used to determine the number states to be retained.
GS_THRESHOLD=1.000000000E-09 !Energy threshold for ground state degeneracy loop up
HWBAND=2.000000000 !half-bandwidth for the bath initialization: flat in -hwband:hwband
LANC_METHOD=arpack !select the lanczos method to be used in the determination of the spectrum. ARPACK (default), LANCZOS (T=0 only), DVDSON (no MPI)
LANC_NSTATES_SECTOR=1 !Initial number of states per sector to be determined.
LANC_NSTATES_TOTAL=1 !Initial number of total states to be determined.
LANC_NSTATES_STEP=2 !Number of states added to the spectrum at each step.
LANC_NCV_FACTOR=10 !Set the size of the block used in Lanczos-Arpack by multiplying the required Neigen (Ncv=lanc_ncv_factor*Neigen+lanc_ncv_add)
LANC_NCV_ADD=0 !Adds up to the size of the block to prevent it to become too small (Ncv=lanc_ncv_factor*Neigen+lanc_ncv_add)
LANC_NITER=512 !Number of Lanczos iteration in spectrum determination.
LANC_NGFITER=300 !Number of Lanczos iteration in GF determination. Number of momenta.
LANC_TOLERANCE=1.000000000E-12 !Tolerance for the Lanczos iterations as used in Arpack and plain lanczos.
LANC_DIM_THRESHOLD=128 !Min dimension threshold to use Lanczos determination of the spectrum rather than Lapack based exact diagonalization.
CG_METHOD=0 !Conjugate-Gradient method: 0=NR, 1=minimize.
CG_GRAD=1 !Gradient evaluation method: 0=analytic (default), 1=numeric.
CG_FTOL=1.000000000E-05 !Conjugate-Gradient tolerance.
CG_STOP=0 !Conjugate-Gradient stopping condition: 0-3, 0=C1.AND.C2, 1=C1, 2=C2 with C1=|F_n-1 -F_n|<tol*(1+F_n), C2=||x_n-1 -x_n||<tol*(1+||x_n||).
CG_NITER=1000 !Max. number of Conjugate-Gradient iterations.
CG_WEIGHT=1 !Conjugate-Gradient weight form: 1=1.0, 2=1/n , 3=1/w_n.
CG_SCHEME=delta !Conjugate-Gradient fit scheme: delta or weiss.
CG_POW=2 !Fit power for the calculation of the Chi distance function as 1/L*|G0 - G0and|**cg_pow
CG_MINIMIZE_VER=F !Flag to pick old/.false. (Krauth) or new/.true. (Lichtenstein) version of the minimize CG routine
CG_MINIMIZE_HH=1.000000000E-04 !Unknown parameter used in the CG minimize procedure.
BATH_TYPE=normal !flag to set bath type: normal (1bath/imp), hybrid(1bath), replica(1replica/imp)
HFILE=hamiltonian !File where to retrieve/store the bath parameters.
IMPHFILE=inputHLOC.in !File read the input local H.
LOGFILE=6 !LOG unit.
19 changes: 7 additions & 12 deletions test/python/hm_bethe.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,10 +47,6 @@ def dens_bethe(x,d):
ed.set_hloc(Hloc)
Nb=ed.get_bath_dimension()
bath = ed.init_solver()
ed.finalize_solver()
print("AAAAA")
ed.init_solver(Nb)
ed.set_hloc(hloc=Hloc)
bath_prev = np.copy(bath)

#DMFT CYCLE
Expand All @@ -71,14 +67,12 @@ def dens_bethe(x,d):
Gmats = np.zeros_like(Smats)
Greal = np.zeros_like(Sreal)

#Compute the local gf:
for i in range(ed.Lmats):
zeta = wm[i]*1j - Smats[0,0,0,0,i]
Gmats[0,0,0,0,i] = sum(Dband/(zeta - Eband))*de

for i in range(ed.Lreal):
zeta = wr[i]+ed.eps*1j - Sreal[0,0,0,0,i]
Greal[0,0,0,0,i] = sum(Dband/(zeta - Eband))*de
zeta = wm*1j - Smats[0,0,0,0,:]
Gmats[0,0,0,0,:] = np.sum(Dband / (zeta[:, np.newaxis] - Eband), axis=1)*de

zeta = wr+ed.eps*1j - Sreal[0,0,0,0,:]
Greal[0,0,0,0,:] = np.sum(Dband / (zeta[:, np.newaxis] - Eband), axis=1)*de

if(rank==0):
np.savetxt('Gloc_iw.dat', np.transpose([wm,Gmats[0,0,0,0,:].imag,Gmats[0,0,0,0,:].real]))
Expand All @@ -95,7 +89,8 @@ def dens_bethe(x,d):
bath = wmixing*bath + (1.0-wmixing)*bath_prev
bath_prev=np.copy(bath)

err,converged=ed.check_convergence(Delta[0,0,0,0,:],ed.dmft_error,1,ed.Nloop)
err,converged=ed.check_convergence(Delta,ed.dmft_error)

ed.finalize_solver()
print("Done...")

0 comments on commit f4be00f

Please sign in to comment.