====== Wrong (missing) ground state Part 2 ====== ;;# asked by [[mailto:riccardo.piombo@gmail.com|Riccardo Piombo]] (2020/01/07 18:57) ;;# == == 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 ~~DISCUSSION|Answers~~