PES

Hi there, it seems I've got a problem in PES spectrum computation: I obtain the correct spectrum of a d10 system if I set all the Coulomb interaction equal to zero. It should be the same also in the interacting case because I'm adding a single hole to the system, but it doesn't happen. Anyone can give me any advice?

Here is the part of the script where I compute the PES Spectrum. If necessary, I can also send you the part of the code in which the Hamiltonian is defined.

A small clarification: the system is a molecule formed by a transition metal and a Ligand atom. I'm considering only the sub-shell 4d of the TM while I'm considering a sub-shell with the same symmetry as the 4d in the Ligand: the latter is formed by the bonding linear combinations of p spin-orbitals.

-- number of fermionic spin-orbitals
NF=20
-- number of bosonic modes
NB=0

-- 4d spin-orbitals in spherical harmonics basis
IndexDn_4d={ 0, 2, 4, 6, 8}
IndexUp_4d={ 1, 3, 5, 7, 9}

-- Ligand-P spin-orbitals in spherical harmonics basis
IndexDn_Ld={10, 12, 14, 16, 18}
IndexUp_Ld={11, 13, 15, 17, 19}

--#########################################  PES Spectrum  #########################################
 
----------- Hamiltonian diagonalization -----------
Npsi_i = 1
psiList_N = Eigensystem(H_tot, StartRestrictions1, Npsi_i)
psiList_N = {psiList_N} 

oppList={H_tot, OppHopLd4d, OppCF_4d, OppCF_Ld, OppUdd, OppULpLp, OppN_4d, OppN_Ld}
print("Done")

print("")
print("---------- Energies in eV ----------")
print("")
print("  #    <E>      <Ehyb>    <CF4d>    <CFLd>     <Udd>    <Upp>     <N_4d>    <N_Ld>");
for i = 1,#psiList_N do
  io.write(string.format("%3i ",i))
  for j = 1,#oppList do
    expectationvalue = Chop(psiList_N[i]*oppList[j]*psiList_N[i])
    io.write(string.format("%9.5f ",expectationvalue))
  end
  io.write("\n")
end

if ToDo['Scalar Product'] then
  print("")
  print("---------- Scalar product <ψ[i]|ψ[j]> ----------")
  print("")
  for i = 1,#psiList_N do
    for j = 1,#psiList_N do
      io.write(string.format("%3i %3i ",i,j))
      scalar_product = Chop(psiList_N[i]*psiList_N[j])
      io.write(string.format("%8.3f ",scalar_product))
    end
    io.write("\n")
  end
end

----------- chemical potential -----------

-- Fermi energy is not defined in the experimental data so we have to fix 
-- it and then freely translate the experimental spectrum along w axis
-- Three possible choices: μ+ , μ- , mean μ

-- μ- : chemical potential to remove a particle
print("")
print("(N-1)-body states energies computation")
StartRestrictions_rem = {NF, NB, {"1111111111 1111111111",nL + (n4d-1),nL + (n4d-1)}}
psiList_Nminus1 = Eigensystem(H_tot, StartRestrictions_rem, 1)
psiList_Nminus1 = {psiList_Nminus1}

E_gs_N = psiList_N[1]*(H_tot)*psiList_N[1]
E_gs_Nminus1 = psiList_Nminus1[1]*(H_tot)*psiList_Nminus1[1]
remove_mu = E_gs_N - E_gs_Nminus1 
print("E_gs_N = ", E_gs_N)
print("E_GS_N-1 = ", E_gs_Nminus1)
print("mu- = ", remove_mu)
print("Done")
print("")

-- μ+ : chemical potential to add a particle
print("(N+1)-body states energies computation")
StartRestrictions_add = {NF, NB, {"1111111111 1111111111",nL + (n4d+1), nL + (n4d+1)}}
psiList_Nplus1 = Eigensystem(H_tot, StartRestrictions_add, 1)
psiList_Nplus1 = {psiList_Nplus1}
-- check for n4d = 10
-- print(psiList_Nplus1) == is nhil
-- E_Gs(N+1) = 0

E_gs_Nplus1 = psiList_Nplus1[1]*(H_tot)*psiList_Nplus1[1]
add_mu = E_gs_Nplus1 - E_gs_N
print("E_gs_N = ", E_gs_N)
print("E_GS_N+1 = ", E_gs_Nplus1)
print("mu+ = ", add_mu)
print("Done")
print("")

-- mean μ   
mean_mu = (add_mu + remove_mu)/2.0
print("mean mu = [(mu+) + (mu-)]/2 = ", mean_mu)
print("")

----------- change-basis matrix -----------

Orbitals = {"4d","Ligand-P_d"}
YtoK = YtoKMatrix(Orbitals)
KtoY = ConjugateTranspose(YtoK)
print("YtoK = ")
print(YtoK)
print("")
print("KtoY =")
print(KtoY)


----------- PES Spectrum Computation -----------

-- Rotate returns M' = R* M R^T
-- the hermitian conjugate to KtoY is YtoK
-- we use rotate to obtain Kubic harmonics

-- We create a list of operators related to the removal and adding spectrum
-- Then the spectral function of each operator in this list is computed
-- the output file contains as many columns as many operators are specified in the list. 
-- The columns can be later combined using a python program to obtain the operators in 
-- the requested representation

-- from T_SW[1] to T_SW[20] there are the operators related to the removal spectrum 
T_SW = {}
for k,v in pairs(IndexUp_4d) do
    T_SW[#T_SW + 1] = Rotate(NewOperator("An", NF, v-1), KtoY) ---> v-1 provides for the down components
    -- here #T_SW increases by one
    T_SW[#T_SW + 1] = Rotate(NewOperator("An", NF, v), KtoY) ---> v provides for the up components
end

print(T_SW)

for k,v in pairs(IndexUp_Ld) do
    T_SW[#T_SW + 1] = Rotate(NewOperator("An", NF, v-1), KtoY)
    T_SW[#T_SW + 1] = Rotate(NewOperator("An", NF, v), KtoY)
end

-- from T_SW[21] to T_SW[40] there are the operators related to the adding spectrum 
for k,v in pairs(IndexUp_4d) do
    T_SW[#T_SW + 1] = Rotate(NewOperator("Cr", NF, v-1), KtoY)
    T_SW[#T_SW + 1] = Rotate(NewOperator("Cr", NF, v), KtoY)
end
for k,v in pairs(IndexUp_Ld) do
    T_SW[#T_SW + 1] = Rotate(NewOperator("Cr", NF, v-1), KtoY)
    T_SW[#T_SW + 1] = Rotate(NewOperator("Cr", NF, v), KtoY)
end

print("Gamma  = ", Gamma)
print("Npoint = ", Npoint)
print("Gamma unb  = ", Gamma_unb)
print("Npoint unb = ", Npoint_unb)

-- CreateSpectra returns G(ω) = <ψ|O*_2 (ω + E0 + iΓ/2 - O_1)^-1 O_2|ψ> 
-- where E0 = <ψ|O_1|ψ>. Here O_2 = Σ_m d_m and O_1 = Htot. 

-- Total spectral wieght
G_SW = CreateSpectra(H_tot, T_SW , psiList_N[1], {{"Emin",E_min}, {"Emax",E_max}, {"NE",Npoint}, {"Gamma",Gamma}} )
--G_SW.Shift(remove_mu) 
G_SW = -(1/math.pi)*G_SW

peaks_SW = CreateSpectra(H_tot, T_SW, psiList_N[1], {{"Emin", E_min}, {"Emax",E_max}, {"NE",Npoint_unb}, {'Gamma',Gamma_unb}} )
--peaks_SW.Shift(remove_mu)
peaks_SW = -(1/math.pi)*peaks_SW

G_SW.Print({{"file","PES_Spectra.dat"}})
peaks_SW.Print({{"file", "Peaks.dat"}})