dofile("../definitions.Quanty")
-- we define an arbitrary operator
Opp = -(2*OppSy + OppLy)*(2*OppSy + OppLy) + (2*OppSy + OppLy) + 0.0000001 * OppLy
-- and a starting density
rho = DensityMatrix(psi1)
-- as well as a temperature needed to average over degenerate states
T = 0.0001
-- calculate the ground-state in mean-field
rhogrd, E, OppMF = MeanFieldGroundState(Opp, rho, T)
-- print the resulting density
print(Chop(rhogrd))
-- print the ground-state energy
print(E)
-- print the Hamiltonian in mean-field, i.e. a potential optimized for the ground-state density
print(Chop(OppMF))
-- lets compare the eigenstates of Opp and OppMF
Npsi=15
-- In order to make sure we have a filling of 2
-- electrons we need to define some restrictions
StartRestrictions = {Nf, Nb, {"111111",2,2}}
-- We now can create the lowest Npsi eigenstates:
psiList = Eigensystem(Opp, StartRestrictions, Npsi)
psiListMF = Eigensystem(OppMF, StartRestrictions, Npsi)
-- We define a list of some operators to look at their expectation values
oppList={Opp, OppMF, OppSy, OppLy}
-- after we've created the eigen states we can look
-- at a list of their expectation values
-- on the left we show the full eigen-states, on the right the eigen-states of the mean-field approximated operator
print(" MF ");
for i=1,#psiList do
for j=1,#oppList do
io.write(string.format("%7.3f ",Chop(psiList[i]*oppList[j]*psiList[i])))
end
io.write(" | ")
for j=1,#oppList do
io.write(string.format("%7.3f ",Chop(psiListMF[i]*oppList[j]*psiListMF[i])))
end
io.write("\n")
end