# 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
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 = {psiList_Nplus1}
print("Done")

E_gs_Nplus1 = psiList_Nplus1[1]*(H_tot)*psiList_Nplus1[1]
print("N plus 1", E_gs_Nplus1)
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
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

, 2020/01/07 20:35

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

psiList_Nminus1 = Eigensystem(H_tot, StartRestrictions_rem, 1)

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

 StartRestrictions_rem = {NF, NB, {"1111111111 1111111111",18 18}}

would work.

Best wishes, Maurits