# PES

asked by Riccardo Piombo (2020/04/30 12:10)

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 = {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]
print("E_gs_N = ", E_gs_N)
print("E_GS_N+1 = ", E_gs_Nplus1)
print("Done")
print("")

-- mean μ
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"}})


, 2020/05/13 07:52

Dear Ricardo,

I have not had the time to look in detail to your files, but if I understood your problem correctly your input file works for the case where you set F0=0, but not for the case where F0 is finite.

Correlations will shift the range of the spectrum. The models we normally use (ligand field theory with the parameterisation in terms of U and Delta do not set the chemical potential to the right value. This is not a problem as we restrict our Hilbert space to a given filling. It does mean that a PES spectrum can end up going from -100 to -80 eV with the Fermi energy / Chemical potential at -100. (Any other arbitrary range of numbers between -1000 and 1000 eV is allowed)

As a first try, what happens if you set the spectrum range from -200 to 200 and the broadening to 10 eV?

PS you can calculate the chemical potential form you model

, 2020/05/13 19:18, 2020/05/13 19:40

Dear Maurits, If there is the possibility to upload some images I will show you the differences, but I do not know how to do it. In particular, the shape of the PES Spectrum changes between the two cases (interacting and non-interacting).

When the interaction is turned on the whole PES Spectrum shifts about 50 eV which is roughly 9*Udd (in my simulation Udd is approximately 5.8 eV). Furthermore, the set of the energy states changes. It can be seen from the position of the peaks found through the diagonalization of the Hamiltonian.

Thanks, Riccardo

, 2020/05/14 09:42

Yes, adding pictures would be nice. (I'm not super happy with docu wiki, especially not with the options to make a forum). If you sent them to me by mail I'll add them to the Forum post.

Maurits

, 2020/05/18 14:20, 2020/05/19 16:17

Sorry for editing completely the message but I found a “strange bug” that could be responsible for the problem. To simplify you can just take a shell d with 10 spin-orbital and this hamiltonian

H = lambda*OppUdd + Bz*(OppLz + 2*OppSz) + zeta_3d*Oppldots + (ed + Udd*(1-lambda))*OppN_d

lambda is a coupling constant to turn on/off the Coulomb repulsion. While diagonalizing a d9 system I obtain random values for the J^2 eigenvalue. Here are some example of the output. Look at the J^2 column.

Run
<F[0]> <F[2]> <F[4]>  <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <J_z>   <E>  | Term
36.000 -1.143 -1.143  0.750  6.000  6.453  0.031  0.716  0.748 274.241 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  7.138 -0.044 -0.022 -0.066 274.238 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  6.625 -0.003  0.570  0.567 274.241 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  8.433 -0.186 -0.894 -1.080 274.232 | ^2D_2.5
36.000 -1.143 -1.143  0.750  6.000  6.467  0.362 -1.226 -0.864 274.241 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  6.902 -0.109  0.857  0.748 274.239 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  6.599 -0.289  0.133 -0.156 274.241 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  5.198  0.068  0.138  0.206 274.248 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  5.288 -0.018 -1.004 -1.022 274.247 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  8.398  0.188  0.730  0.919 274.232 | ^2D_2.5

Run
<F[0]> <F[2]> <F[4]>  <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <J_z>   <E>  | Term
36.000 -1.143 -1.143  0.750  6.000  5.319 -0.043 -0.137 -0.179 274.247 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  8.269  0.221  1.206  1.427 274.232 | ^2D_2.5
36.000 -1.143 -1.143  0.750  6.000  6.781 -0.251 -0.923 -1.174 274.240 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  7.807  0.295  0.964  1.259 274.235 | ^2D_2.5
36.000 -1.143 -1.143  0.750  6.000  6.223 -0.177  0.198  0.021 274.243 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  7.982 -0.161 -0.468 -0.628 274.234 | ^2D_2.5
36.000 -1.143 -1.143  0.750  6.000  5.062 -0.062  0.422  0.359 274.248 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  7.098 -0.120 -1.011 -1.131 274.238 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  6.761 -0.085 -0.500 -0.585 274.240 | ^2D_2
36.000 -1.143 -1.143  0.750  6.000  6.199  0.382  0.249  0.631 274.243 | ^2D_2

I could have shown you countless others. At each run that column changes randomly BUT If I repeat the same calculation for a d^1 system (equal to the d^9 one thanks to particle-hole symmetry) I instead obtain the right multiplets:

<F[0]> <F[2]> <F[4]>  <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <J_z>   <E>  | Term
0.000  0.000  0.000  0.750  6.000  3.750  0.300 -1.800 -1.500 -0.015 | ^2D_1.5
0.000  0.000  0.000  0.750  6.000  3.750  0.100 -0.600 -0.500 -0.015 | ^2D_1.5
0.000  0.000  0.000  0.750  6.000  3.750 -0.100  0.600  0.500 -0.015 | ^2D_1.5
0.000  0.000  0.000  0.750  6.000  3.750 -0.300  1.800  1.500 -0.015 | ^2D_1.5
0.000  0.000  0.000  0.750  6.000  8.750 -0.500 -2.000 -2.500  0.010 | ^2D_2.5
0.000  0.000  0.000  0.750  6.000  8.750 -0.300 -1.200 -1.500  0.010 | ^2D_2.5
0.000  0.000  0.000  0.750  6.000  8.750 -0.100 -0.400 -0.500  0.010 | ^2D_2.5
0.000  0.000  0.000  0.750  6.000  8.750  0.100  0.400  0.500  0.010 | ^2D_2.5
0.000  0.000  0.000  0.750  6.000  8.750  0.300  1.200  1.500  0.010 | ^2D_2.5
0.000  0.000  0.000  0.750  6.000  8.750  0.500  2.000  2.500  0.010 | ^2D_2.5

It's kinda weird, don't you think? It involves the d9 system only!!

, 2020/05/19 16:18, 2020/05/19 16:19

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"}})
, 2020/05/25 17:29

I made some tests and I've found that while diagonalizing a d10 system if zeta spin-orbit is too small with respect to F0 then states with different J mixes and j is not still a good quantum number. The problem is resolved by putting a larger spin-orbit but if my system has a spin-orbit coupling of the order of 0.01 eV and an F0 of the order of 6 eV for the 4d levels and Quanty is not able to diagonalize the Hamiltonian correctly with these values, how can I do it? Is there any variable related to the number of iterations involved in the diagonalization?

, 2020/05/25 22:04

Dear Riccardo,

Quanty stops the calculation when a given accuracy is reached. Determining if this is the case is not always easy. The difference between 10.001 and 10.002 might seem small, but is the same as the difference between 0.001 and 0.002 where it is a factor of two. For your d9 case the energy difference between the two SOC states is on the sub per mile level. This is due to a shift of the energy of all your states. The code should recognise this, but in your case did not. I will have a look to see if I can improve this.

At the same time you can tell the code the accuracy of the calculation you want

Change line 117 from

psi_i = Eigensystem(H, Restrictions_i, Nstates_i)    --> diagonalization of H

to

psi_i = Eigensystem(H, Restrictions_i, Nstates_i,{{"zero",1E-10}})    --> diagonalization of H

and the problem should disappear.

Sorry for that

Maurits

, 2020/05/13 07:53, 2020/05/13 07:57
-- CreateSpectra returns G(ω) = <ψ|O*_2 (ω + E0 + iΓ/2 - O_1)^-1 O_2|ψ>

I've seen you noticed that Quanty / Lua can handle more than ASCII characters….. (You can now generate the most obfuscated code by using the different space characters in unicode as the different variables. On top of that redefine all functions by similar variables made of space characters and you can have an seemingly empty file doing your calculations…)

, 2020/05/25 22:13, 2020/05/25 22:14

PS If you would leave the verbosity to its standard value you would get this extra output With standard value for convergence

=========================================================
N-body states energy computation
=========================================================
Start of BlockGroundState. Converge 10 states to an energy with relative variance smaller than  1.490116119384766E-06

Start of BlockOperatorPsiSerialRestricted
Start of BlockGroundState. Converge 10 states to an energy with relative variance smaller than  1.490116119384766E-06

Start of BlockOperatorPsiSerial

With better convergence

=========================================================
N-body states energy computation
=========================================================
Start of BlockGroundState. Converge 10 states to an energy with relative variance smaller than  1.000000000000000E-10

Start of BlockOperatorPsiSerialRestricted
Outer loop   1, Number of Determinants:        10        10 last variance  1.559885131428018E-04
Start of BlockOperatorPsiSerialRestricted
Start of BlockGroundState. Converge 10 states to an energy with relative variance smaller than  1.000000000000000E-10

Start of BlockOperatorPsiSerial

Quanty decided that the random starting vectors were converged enough for the accuracy asked for and never started to diagonalise… When I start to diagonalise I apply a shift, and after the shift the variance is no longer below the threshold.