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