PES

Here is the script. You can run it to check this strange thing

Input.in is a file with l and Nelec

Verbosity(0)

-- define function which will print the Term symbol |2S+1 L_J> according to L^2 S^2 and J^2
function LSJ2term (S2,L2,J2)
  L = floor(0.5 * (sqrt(abs(L2) * 4 +1) - 1) + 0.5)
  if     L == 0   then str1="S"
  elseif L == 1   then str1="P"
  elseif L == 2   then str1="D"
  elseif L == 3   then str1="F"
  elseif L == 4   then str1="G"
  elseif L == 5   then str1="H"
  elseif L == 6   then str1="I"
  elseif L == 7   then str1="K"
  elseif L == 8   then str1="L"
  elseif L == 9   then str1="M"
  elseif L == 10  then str1="N"
  elseif L == 11  then str1="O"
  elseif L == 12  then str1="Q"
  elseif L == 13  then str1="R"
  elseif L == 14  then str1="T"
  elseif L == 15  then str1="U"
  elseif L == 16  then str1="V"
  elseif L == 17  then str1="W"
  elseif L == 18  then str1="X"
  elseif L == 19  then str1="Y"
  elseif L == 20  then str1="Z"
  else 
    str1="?"
  end
  mult = floor(sqrt(S2 * 4 +1)+0.5)
  twoJ = floor((sqrt(abs(J2) * 4 +1) - 1) + 0.5)/2
  --print(twoJ)
  return "^"..mult..str1.."_"..twoJ
end

print("")
print("=================================================================================")
print("--> This program computes the atomic VPES spectrum. Don't forget to set the chosen") 
print("--> config in input.in")
print("=================================================================================")
print("")

dofile("input.in")

NF = 2*(2*l+1)   --> number of spin-orbitals
NB = 0           --> number of bosonic modes

Nk = l+1   --> multipole expansion order

-- spin-orbitals in |lm> basis
IndexDn={}
IndexUp={}
 
j=1
for i=0,NF-1,2 do
  IndexDn[j]=i
  IndexUp[j]=i+1
  j=j+1
end

-- define few Operator on this basis
OppSz   = NewOperator("Sz"   , NF, IndexUp, IndexDn)
OppLz   = NewOperator("Lz"   , NF, IndexUp, IndexDn)
OppSsqr = NewOperator("Ssqr" , NF, IndexUp, IndexDn)
OppLsqr = NewOperator("Lsqr" , NF, IndexUp, IndexDn)
OppJz   = NewOperator("Jz"   , NF, IndexUp, IndexDn)
OppJsqr = NewOperator("Jsqr" , NF, IndexUp, IndexDn)
Oppldots= NewOperator("ldots", NF, IndexUp, IndexDn)

A = 7.4
B = 0.15
C = 0.58

F0  = A + 7.0*C/5.0
F2  = 49.0*B + 7*C
F4  = 441.0*C/35.0
Udd = F0 - (F2 + F4)*2.0/63.0 
print("Udd = ", Udd)
print("F[0] = ", F0)
print("F[2] = ", F2)
print("F[4] = ", F4)

OppFk={}
OppFk[1] = NewOperator("U", NF, IndexUp, IndexDn, {1,0,0})
OppFk[2] = NewOperator("U", NF, IndexUp, IndexDn, {0,1,0})
OppFk[3] = NewOperator("U", NF, IndexUp, IndexDn, {0,0,1})

OppUdd = F0*OppFk[1] --+ F2*OppFk[2] + F4*OppFk[3]
OppUdd.Name = "d-d Coulomb repulsion"

OppNUp_d = NewOperator("Number", NF, IndexUp, IndexUp, {1,1,1,1,1})
OppNDn_d = NewOperator("Number", NF, IndexDn, IndexDn, {1,1,1,1,1})
OppN_d = OppNUp_d + OppNDn_d
OppN_d.Name = "d Number operator"

zeta_3d = 0.01
Bz = 0.000001*EnergyUnits.Tesla.value
lambda = 1
ed=0.

-- when lambda = 0 and ed = 0, <H> is approximately N_d*Udd
H = lambda*OppUdd + Bz*(OppLz + 2*OppSz) + zeta_3d*Oppldots + (ed + Udd*(1-lambda))*OppN_d 
    

print("")
print("=========================================================")
print("N-body states energy computation")
print("=========================================================")

Restrictions_i = {NF, NB, {"1111111111", Nelec, Nelec}}  --> inital restrictions for diagonalization

----------- N-Hamiltonian diagonalization -----------
Nstates_i = Binomial(NF,Nelec)                       --> number of inital states
psi_i = Eigensystem(H, Restrictions_i, Nstates_i)    --> diagonalization of H
--psi_i = {psi_i}

E_gs_i = psi_i[1]*(H)*psi_i[1]  --> gs energy of the N-body system

print("")
print("--> After N-body H diagonalization")
print("Nelec = ", Nelec)
print("")

for k=1,Nk do
  io.write(string.format("<F[%i]> ",2*(k-1)))
end

oppList = {OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz, OppJz, H}

print(" <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <J_z>   <E>  | Term");
Res={}
ResF={}
for i=1,Nstates_i do
  for k=1,Nk do
    ResF[k] = psi_i[i] * OppFk[k] * psi_i[i]
    io.write(string.format("%6.3f ",ResF[k]))
  end
  for j=1,7 do  --> 7 is the number of entries in OppList
    Res[j] = psi_i[i] * oppList[j] * psi_i[i]
    io.write(string.format("%6.3f ",Res[j]))
  end
  io.write(string.format("| %s",LSJ2term(Res[1],Res[2],Res[3])))
  io.write("\n")
end

print("")
print("=========================================================")
print("(N-1)-body states energy computation")
print("=========================================================")

Restrictions_f = {NF, NB, {"1111111111", Nelec-1, Nelec-1}}

----------- (N-1)-Hamiltonian diagonalization -----------
Nstates_f = Binomial(NF,Nelec-1)
print(Nstates)
psi_f = Eigensystem(H, Restrictions_f, Nstates_f)

E_gs_f = psi_f[1]*(H)*psi_f[1]  --> gs energy of the (N-1)-body system
DeltaE_gs = E_gs_i - E_gs_f

print("")
print("--> After (N-1)-body H diagonalization")
print("Nelec = ", Nelec-1)
print("")

for k=1,Nk do
  io.write(string.format("<F[%i]> ",2*(k-1)))
end

oppList = {OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz, OppJz, H}

print(" <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <J_z>   <E>  | Term");
Res_2={}
ResF_2={}
for i=1,Nstates_f do
  for k=1,Nk do
    ResF_2[k] = psi_f[i] * OppFk[k] * psi_f[i]
    io.write(string.format("%6.3f ",ResF_2[k]))
  end
  for j=1,7 do  --> 7 is the number of entries in OppList
    Res_2[j] = psi_f[i] * oppList[j] * psi_f[i]
    io.write(string.format("%6.3f ",Res_2[j]))
  end
  io.write(string.format("| %s",LSJ2term(Res_2[1],Res_2[2],Res_2[3])))
  io.write("\n")
end

print("")
print("E_gs_i = ", E_gs_i)
print("E_gs_f = ", E_gs_f)
print("Delta E_gs = ", DeltaE_gs)

Orbitals = {"3d"}
YtoK = YtoKMatrix(Orbitals)
KtoY = ConjugateTranspose(YtoK)

T_PES = {}

-- from T_PES[1] to T_PES[20] there are the operators related to the removal spectrum 
-- Here is the TM removal part
for k,v in pairs(IndexUp) do
  T_PES[#T_PES + 1] = Rotate(NewOperator("An", NF, v-1), KtoY) ---> v-1 provides for the down components
  T_PES[#T_PES + 1] = Rotate(NewOperator("An", NF, v), KtoY) ---> v provides for the up components
end
--print(T_PES)

E_max =  100
E_min = -100
Gamma =  1.6
Gamma_unb = 0.01
Npoint = 100*(E_max-E_min)/Gamma
Npoint_unb = 10*(E_max-E_min)/Gamma_unb

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

-- Total spectral wieght
PESSpectra = CreateSpectra(H, T_PES , psi_i[1], {{"Emin",E_min}, {"Emax",E_max}, {"NE",Npoint}, {"Gamma",Gamma}} )
--PESSpectra.Shift(remove_mu) 
PESSpectra = -(1/math.pi)*PESSpectra

peaks_PESSpectra = CreateSpectra(H, T_PES, psi_i[1], {{"Emin", E_min}, {"Emax",E_max}, {"NE",Npoint_unb}, {'Gamma',Gamma_unb}} )
--peaks_PESSpectra.Shift(remove_mu)
peaks_PESSpectra = -(1/math.pi)*peaks_PESSpectra

PESSpectra.Print({{"file","PES_Spectra.dat"}})
peaks_PESSpectra.Print({{"file", "Peaks.dat"}})