add test script

test/python/hamiltonian.restart

@@ -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

test/python/hm_bethe.py

@@ -0,0 +1,97 @@
+import numpy as np
+import scipy as sp
+from edipy2 import global_env as ed
+import mpi4py
+from mpi4py import MPI
+import os,sys
+
+#INIT MPI 
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+print("I am process",rank,"of",comm.Get_size())
+master = (rank==0)
+
+def dens_bethe(x,d):
+  root=(1-(x/d)**2)+0.j
+  root=np.sqrt(root)
+  dens_bethe=(2/(np.pi*d))*root
+  return dens_bethe.real
+
+#READ ED INPUT:
+ed.read_input("inputED.conf")
+if(ed.Nspin!=1 or ed.Norb!=1):
+    print("This test code is for Nspin=1 + Norb=1.")
+    ed.Nspin=1
+    ed.Norb=1
+Le      = 1000
+wmixing = 0.3
+wband   = 1.0
+
+#BUILD Density of States:
+Eband,de = np.linspace(-wband,wband,Le,retstep=True)
+Dband = dens_bethe(Eband,wband)
+
+#BUILD frequency array:
+wm = np.pi/ed.beta*(2*np.arange(ed.Lmats)+1)
+wr = np.linspace(ed.wini,ed.wfin,ed.Lreal)
+
+
+#Allocate minimal required arrays:
+#ALL functions must have shape [Nspin,Nspin,Norb,Norb(,L)]:
+Smats=np.zeros((ed.Nspin,ed.Nspin,ed.Norb,ed.Norb,ed.Lmats),dtype='complex',order='F')
+Gmats=np.zeros((ed.Nspin,ed.Nspin,ed.Norb,ed.Norb,ed.Lmats),dtype='complex',order='F')
+Sreal=np.zeros((ed.Nspin,ed.Nspin,ed.Norb,ed.Norb,ed.Lreal),dtype='complex',order='F')
+Greal=np.zeros((ed.Nspin,ed.Nspin,ed.Norb,ed.Norb,ed.Lreal),dtype='complex',order='F')
+Delta=np.zeros((ed.Nspin,ed.Nspin,ed.Norb,ed.Norb,ed.Lmats),dtype='complex',order='F')
+Hloc =np.zeros((ed.Nspin,ed.Nspin,ed.Norb,ed.Norb),dtype='float',order='F')
+
+
+
+#SETUP SOLVER
+Nb=ed.get_bath_dimension()
+bath = np.zeros(Nb,dtype='float',order='F')
+ed.init_solver(bath)
+bath_prev = np.copy(bath)
+
+#DMFT CYCLE
+converged=False;iloop=0
+while (not converged and iloop<ed.Nloop ):
+    iloop=iloop+1
+    print("DMFT-loop:",iloop,"/",ed.Nloop)
+
+    #Solve the EFFECTIVE IMPURITY PROBLEM (first w/ a guess for the bath)
+    ed.solve(bath,Hloc)
+
+
+    #Get Self-energies
+    ed.get_sigma_matsubara(Smats)
+    ed.get_sigma_realaxis(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
+
+    if(rank==0):
+        np.savetxt('Gloc_iw.dat', np.transpose([wm,Gmats[0,0,0,0,:].imag,Gmats[0,0,0,0,:].real]))
+        np.savetxt('Gloc_realw.dat', np.transpose([wr,Greal[0,0,0,0,:].imag,Greal[0,0,0,0,:].real]))
+
+
+    #Get the Delta function and FIT:    
+    Delta[0,0,0,0,:] = 0.25*wband*Gmats[0,0,0,0,:]
+    if(rank==0):
+        np.savetxt('Delta_iw.dat', np.transpose([wm,Delta[0,0,0,0,:].imag,Delta[0,0,0,0,:].real]))
+    ed.chi2_fitgf(Delta,bath,ispin=1,iorb=1)
+
+    if(iloop>1):
+        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)
+
+print("Done...")
+

test/python/inputED.conf

@@ -0,0 +1,78 @@
+ LE=1000                                        !
+ WBETHE=1.000000000,1.000000000,1.000000000,1.000000000,1.000000000 !
+ DBETHE=0.d0,0.d0,0.d0,0.d0,0.d0               !
+ WMIXING=5.000000000E-01                       !
+ 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=2                                  !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.