STRANGER THINGS

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 by 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 ciclycally three different value for the GS energy which differ from one another by 1.6 eV 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 [1] B = 0.09 C = 0.54

-- Slater-Koster integrals for monopole and -- multipole part of Coulomb interaction [2] 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 = 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_2_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 print(”“) print(”(N-1)-body states energies computation“) StartRestrictions_rem = {NF, NB, {“1111111111 0000000000”,n4d-1,n4d-1}, {“0000000000 1111111111”,nL,nL}} 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”,n4d+1,n4d+1}, {“0000000000 1111111111”,nL,nL}} 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) ”