# Wrong (missing) ground state Part 2

asked by 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

## Answers

Dear Riccardo

I'm not able to reproduce the problem, (we did implement better convergence tests lately) but guess I know where the problem comes from.

The question is related to a similar problem for finding the ground-state of a $d^6$ configuration, so I changed the title to refer to this question as well.

We use iterative methods to find the ground-state. If the state you are in at some point during the iteration has no overlap with the true ground-state you will not find the ground-state, but the lowest state of the same irreducible representation as the state you are in at the moment.

For your system with two holes, the ground-state is not formally $d^8$, but $d^9L^9$. I.e. you have a charge transfer system and remove one electron from the Ligands.

In the first step the line

will converge the system to the ground-state with the start restrictions on. I.e. the ground-state of a $d^8$ ion ($^3F$). As your crystal-field is zero this is a 21 fold degenerate state and your values of $\vec{L}$ and $\vec{S}$ will be some random number that differs each run you take.

In the next step you allow hopping, but already pre-selected a given irreducible representation and will not be able to reach all.

The solution, use less restrictive starting restrictions. Using

would work.

Best wishes, Maurits