Wrong (missing) ground state Part 2

Hi every one, I wrote a simple code to calculate the chemical potential for the addition and for the removal of a particle. The system under consideration is a planar D4h structure composed of a transition metal surrounded by 4 non-metal ligands. The problem is the following: when computing the ground state of the N-1 system Quanty seems to output cyclically three different value for the GS energy which differs by almost 1.6 eV in the extreme case while for the N and N+1 case it doesn't change (as it must be!). Here is the code could someone copy and paste it on his pc and just run it too check this strange thing!?

Verbosity(0)

NF=20
NB=0

-- Silver 4d spin-orbitals
IndexDn_4d={ 0, 2, 4, 6, 8}
IndexUp_4d={ 1, 3, 5, 7, 9}
-- Fluorine 2p spin-orbitals
IndexDn_Ld={10,12,14,16,18}
IndexUp_Ld={11,13,15,17,19}

-- number of electrons in d-shell and Ligand-P shell
n4d = 9
nL = 10

Udd = 5.0
Delta_pd =  1.5
Tpp =  0.4      
TdLb1g = 2.6    

-- Racah parameters 
B = 0.09
C = 0.54

-- Slater-Koster integrals for monopole and 
-- multipole part of Coulomb interaction 
F2dd = 49.0*B + 7*C
F4dd = 441.0*C/35.0
F0dd = Udd + (F2dd + F4dd)*2.0/63.0  

-- Onsite splitting of the 4d shell 
Eda1g =  0.0
Edb1g =  0.0
Edb2g =  0.0
Edeg  =  0.0

-- Onsite splitting of the Ligand P shell 
ELa1g = -Delta_pd - Tpp
ELb1g = -Delta_pd + Tpp
ELb2g = -Delta_pd - Tpp
ELeg  = -Delta_pd

-- Hopping between the Ligand and 4d shell 
TdLa1g = TdLb1g/sqrt(3.0)
TdLb2g = TdLb1g/2.0
TdLeg  = TdLb1g/(2.0*sqrt(2.0))

-- turning U and Delta to onsite energies 
Delta_ZSA = Delta_pd + Tpp/5.0
eL = (Udd*n4d*(1-n4d)/2.0 - n4d*(Delta_ZSA -Udd*n4d))/(n4d+nL)
ed = eL + Delta_ZSA - Udd*n4d 

-- D4h Crystal field on 4d states
Akm_4d = {{0, 0, ed} ,
{2, 0, Eda1g + (-1)*(Edb1g) + (-1)*(Edb2g) + Edeg} ,
{4, 0, (3/10)*((6)*(Eda1g) + Edb1g + Edb2g + (-8)*(Edeg))} ,
{4,-4, (3/2)*((sqrt(7/10))*(Edb1g + (-1)*(Edb2g)))} ,
{4, 4, (3/2)*((sqrt(7/10))*(Edb1g + (-1)*(Edb2g)))} }

OppCF_4d = NewOperator("CF", NF, IndexUp_4d, IndexDn_4d, Akm_4d)
OppCF_4d.Name = "Crystal Field on 4d-states"

-- D4h Crystal field on Ligand-P states
Akm_L = {{0, 0, eL},
{2, 0, ELa1g + (-1)*(ELb1g) + (-1)*(ELb2g) + ELeg} ,
{4, 0, (3/10)*((6)*(ELa1g) + ELb1g + ELb2g + (-8)*(ELeg))} ,
{4,-4, (3/2)*((sqrt(7/10))*(ELb1g + (-1)*(ELb2g)))} ,
{4, 4, (3/2)*((sqrt(7/10))*(ELb1g + (-1)*(ELb2g)))} }

OppCF_Ld = NewOperator("CF", NF, IndexUp_Ld, IndexDn_Ld, Akm_L)
OppCF_Ld.Name = "Crystal Field on Ligand-P states"

-- Hopping t_pd
Akm_dL = {{0, 0, (1/5)*(TdLa1g + TdLb1g + TdLb2g + (2)*(TdLeg))} ,
         {2, 0, TdLa1g + (-1)*(TdLb1g) + (-1)*(TdLb2g) + TdLeg} ,
         {4, 0, (3/10)*((6)*(TdLa1g) + TdLb1g + TdLb2g + (-8)*(TdLeg))} ,
         {4,-4, (3/2)*((sqrt(7/10))*(TdLb1g + (-1)*(TdLb2g)))} ,
         {4, 4, (3/2)*((sqrt(7/10))*(TdLb1g + (-1)*(TdLb2g)))} }


OppHopLd4d  = NewOperator("CF", NF, IndexUp_4d,IndexDn_4d, IndexUp_Ld,IndexDn_Ld,Akm_dL) +  
              NewOperator("CF", NF, IndexUp_Ld,IndexDn_Ld, IndexUp_4d,IndexDn_4d,Akm_dL)
OppHopLd4d.Name = "P-d hybridization"

-- Udd repulsion 
OppF0 = NewOperator("U", NF, IndexUp_4d, IndexDn_4d, {1,0,0})
OppF2 = NewOperator("U", NF, IndexUp_4d, IndexDn_4d, {0,1,0})
OppF4 = NewOperator("U", NF, IndexUp_4d, IndexDn_4d, {0,0,1})
OppUdd = F0dd*OppF0 + F2dd*OppF2 + F4dd*OppF4

-- Silver hamiltonian
H_Ag = OppUdd + OppCF_4d 

-- Fluorine hamiltonian
H_F = OppCF_Ld

-- Ag + F Hamiltonian
H_Cluster = H_Ag + H_F
H_tot = H_Cluster + (-1)*OppHopLd4d

-- Restrictions = {filling of nd electrons in 4d states, filling of nL electrons in Ligand-P states}
print("")
print("N-body states energies computation")
Npsi_i = 1
StartRestrictions1 = {NF, NB, {"1111111111 0000000000",n4d,n4d}, {"0000000000 1111111111",nL,nL}}
psiList_N = Eigensystem(H_tot, StartRestrictions1, Npsi_i)
psiList_N = {psiList_N}
print("Done")

E_gs_N = psiList_N[1]*(H_tot)*psiList_N[1]
print("")
print("Egs N particle = ", E_gs_N)


-- initialization of the variables
remove_mu = 0
add_mu  = 0
mean_mu = 0
E_gs_Nminus1 = 0
E_gs_Nplus1 =  0

-- chemical potential μ- to remove a particle
-- ###########################################################
-- THE PROBLEM SEEMS TO BE HERE
print("")
print("(N-1)-body states energies computation")
StartRestrictions_rem = {NF, NB, {"1111111111 0000000000",8,8}, {"0000000000 1111111111",10,10}}
psiList_Nminus1 = Eigensystem(H_tot, StartRestrictions_rem, 1)
psiList_Nminus1 = {psiList_Nminus1}
print("Done")

E_gs_Nminus1 = psiList_Nminus1[1]*(H_tot)*psiList_Nminus1[1]
remove_mu = E_gs_N - E_gs_Nminus1  
print("N minus 1", E_gs_Nminus1)
print("μ- = ", remove_mu)
print("")
-- ###########################################################

-- chemical potential μ+ to add a particle
print("(N+1)-body states energies computation")
StartRestrictions_add = {NF, NB, {"1111111111 0000000000",10,10}, {"0000000000 1111111111",10,10}}
psiList_Nplus1 = Eigensystem(H_tot, StartRestrictions_add, 1)
psiList_Nplus1 = {psiList_Nplus1}
print("Done")

E_gs_Nplus1 = psiList_Nplus1[1]*(H_tot)*psiList_Nplus1[1]
add_mu = E_gs_Nplus1 - E_gs_N
print("N plus 1", E_gs_Nplus1)
print("μ+ = ", add_mu)
print("")

-- mean μ  --> Fermi energy is not defined in the experimental values
-- so we have to fix it and then freely translate the experimental spectrum along w axis
mean_mu = (add_mu + remove_mu)/2.0
print("mean μ = [(μ+) + (μ-)]/2 = ", mean_mu)

P.S: after some simulations, the problem seems to appear only when n4d-1 is equal to 8 the GS energy becomes -15.628822671742, -14.433720341628 or -14.042673360744