Mean field theory for fermions: Hartree Fock approximation

DoDebug=True;

Take[{a,b,c},4]

Define a function to do one cycle of the HF calculation.  Input parameters are the number of electrons of each spin (assumed equal),  U, and the densities for opposite spins.

Cycle[nk_,U_,nn_]:=Block[{n,i,vals,vecs,ord,ham},
n=Length[nn];(* the length of the chain *)
ham=Array[0&,{n,n}];
For[i=1,i<n,i++,
ham[[i,i+1]]=-1.;
ham[[i+1,i]]=-1.; (* hopping*)
];
For[i=1,i<=n,i++,
ham[[i,i]]=U nn[[i]]; (* MF potential energy *)
];
{energ,vecs}=Eigensystem[N[ham]];
ord=Ordering[energ]; 
energ=Take[energ[[ord]],nk]; (* Sorted energies up to nk *)
vecs=Take[vecs[[ord]],nk];     (* Sorted vectors up to nk *)
If[DoDebug,
Print[energ];
Print[vecs.ham.Transpose[vecs]//Chop//MatrixForm ]; 
];
{Apply[Plus,energ],
Apply[Plus,vecs*Conjugate[vecs]]}//Chop
(* return energy and new average densities *)
]

vec=Array[0&,20];
DoDebug=False;
For[i=1,i<55,i++,
(* Print[vec];*)
If[i>50,DoDebug=True];
{energ,vec}= Cycle[5,5,vec];
Print[energ];]

ListPlot[vec]

{energ,vec}= Cycle[5,5,vec];

ListPlot[vec]

This is the actual HF loop (last argument is the lattice size)

DoHF[nk_,U_,n_]:=Block[{nn,i,eup,edn,nup,ndn,eold,ee},
nup=Array[0&,n];(* initial density *)
nup[[1]]=nk;
ee=100;
For[i=1,i<500,i++,
{edn,ndn}= Cycle[nk,U,nup];
{eup,nup}=Cycle[nk,U,ndn];
eold=ee;
ee=eup+edn-U nup.ndn;
If[Mod[i,10]==0, Print["U=",U,"  E=",ee]]; (* print out energy once in a while *)
If[Abs[ee-eold]<0.00001,Break[]];
];
{ee,nup,ndn}
]

DoDebug=False;

DoHF[10,1.0,20]

0.49914780659882435`

ListPlot[vec]

<<"ErrorBarPlots`"

For[UU=0.0,UU<10,UU+=1,{energ,nup,ndn}=DoHF[10,UU,80];Print[ListPlot[{nup,ndn,nup+ndn},PlotLabel->"U="<>ToString[UU]<>" E="<>ToString[energ],PlotStyle->{{Hue[0.3]},{Hue[0.5]},{Hue[0.8]}},PlotRange->All,Joined->True]]]

UU=1.0;
For[nel=1,nel≤80,nel++,{energ,nup,ndn}=DoHF[nel,UU,80];Print[ListPlot[{nup,ndn,nup+ndn},PlotLabel->"U="<>ToString[UU]<>"   \!\(\*SubscriptBox[\(n\), \(el\)]\)="<>ToString[nel]<>" E="<>ToString[energ/(2 nel)],PlotStyle->{{Hue[0.3`]},{Hue[0.5`]},{Hue[0.8`]}},PlotRange->All,Joined->True]]]

Compare energies with the Exact Diagonalization

UU=2.2;{energ,nup,ndn}=DoHF[4,UU,5]

HH=Hubbard[8,0,5,UU];

HH//Eigenvalues//Sort

(Eigenvalues[HH,5,Method->{Arnoldi,Shift->-8.0000}]//Sort)

DoDebug=False;

lis={};nelup=3;
For[i=0,i<8.,i+=1,
UU=i;Eh=DoHF[nelup,UU,6][[1]];
HH=Hubbard[2 nelup,0,6,UU]; 
Ee=(Eigenvalues[HH,5,Method->{Arnoldi,Shift->-8.0000}]//Sort)[[1]];
AppendTo[lis,{i,Eh,Ee}]];lis

ListPlot[{Transpose[{Transpose[lis][[1]],Transpose[lis][[2]]}],Transpose[{Transpose[lis][[1]],Transpose[lis][[3]]}]}]

ListPlot[Transpose[{Transpose[lis][[1]],Transpose[lis][[2]]-Transpose[lis][[3]]}],PlotRange->All]

HH=Hubbard[6,0,6,UU];