Hubbard model Hamiltonian

Use functions BitAnd,  BitShiftRight, BitShiftLeft

PP[1313]

BitShiftRight[1313,2]

PP[BitShiftRight[1313,2]]

Define the state as base-4 number.  Convenient functions

Pretty print the state (length 5 by default)

PP[v_,L_:5]:=BaseForm[PaddedForm[v,L,NumberPadding->{"0","0"}],4]

PP[1313,10]

Convenient on-site testing functions

IsSinglet[v_,n_]:=(BitAnd[BitShiftRight[v,2n],3]==3)

Check this function

For[v=32,v<64,v++,Print[PP[v]," ",IsSinglet[v,2]]]

IsUp[v_,n_]:=(BitAnd[BitShiftRight[v,2n],3]==2)

IsDn[v_,n_]:=(BitAnd[BitShiftRight[v,2n],3]==1)

IsHole[v_,n_]:=(BitAnd[BitShiftRight[v,2n],3]==0)

Exercise all functions

For[i=0,i<64,i++,Print["i=",PP[i],"  ",IsSinglet[i,1],"  ",IsUp[i,1],"  ",IsDn[i,1],"  ",IsHole[i,1]]]

Electron Hopping with closed boundary conditions

This function returns the new state vector or (-1) if the move fails.

HopRight[v_,n_,spin_]:=Block[{j},(* spin must be 2 for Up and 1 for Dn *)
If[n==0,Return[-1]]; (* no sites to the right *)
If[BitAnd[BitShiftRight[v,2n-2],spin]≠0,Return[-1]];  (* there is already spin there *)
If[BitAnd[BitShiftRight[v,2n],spin]≠spin,Return[-1]];                     (* no spin at this location *)
       v-BitShiftLeft[spin,2n]+BitShiftLeft[spin,2n-2]
]

Check the function: try hopping spin-down (spin=1) to the right from site 2

For[i=0,i<31,i++,i1=HopRight[i,2,1]; Print[PP[i],"  ",If[i1≠-1,PP[i1],-1]]]

HopLeft[v_,n_,spin_]:=Block[{},
If[BitAnd[BitShiftRight[v,2n+2],spin]≠0,Return[-1]];  (* there is already spin there *)
If[BitAnd[BitShiftRight[v,2n],spin]≠spin,Return[-1]];                     (* no spin at this location *)
        v-BitShiftLeft[spin,2n]+BitShiftLeft[spin,2n+2]
	(* new state=remove spin at position n and add spin at position n+1 *)
]

Check the function

For[i=0,i<31,i++,i1=HopLeft[i,1,1]; Print[PP[i],"  ",If[i1≠-1,PP[i1],-1]]]

This function helps with computing the Fermi factor.  Count total number of electrons to the right of position n (exclusive).

Fermi factor correction

Fermi factor for spin down moving to the right, assuming the move is allowed.  (-1) only if there is already spin up at the new location.

FermiFactorDnRight[v_,n_]:=(* spin up on the right *)
If[BitAnd[BitShiftRight[v,2n-2],2]==2,-1,1]

Fermi factor for spin up moving to the right, assuming allowed move.  (-1) if there is a spin down at the old location:

FermiFactorUpRight[v_,n_]:= (* spin down at this location *)
If[BitAnd[BitShiftRight[v,2n],1]==1,-1,1]

FermiFactorUpLeft[v_,n_]:= (* spin down on the left *)
If[BitAnd[BitShiftRight[v,2n+2],1]==1,-1,1]

FermiFactorDnLeft[v_,n_]:= (* spin up on this location *)
If[BitAnd[BitShiftRight[v,2n],2]==2,-1,1]

For[i=0,i<63,i++,i1=HopLeft[i,1,1]; Print[PP[i],"  ",If[i1≠-1,PP[i1],-1],"  ",If[i1≠-1,FermiFactorDnLeft[i,1]," "]]]

Potential Energy: Compute the number of singlets in the chain of length n

{PP[15*4],DigitCount[15*4,4,{3}]}

PP[BitShiftLeft[1,4]-1]

CountSinglets[v_,n_]:=DigitCount[BitAnd[v,BitShiftLeft[1,2n]-1],4,{3}][[1]]

For[i=0,i<31,i++,Print[PP[i],"  ",CountSinglets[i,4]]];

We can now write the function for creating the Hubbard Hamiltonian.  To make it more efficient, one needs to separate the matrix into independent blocks conserving (a) number of particles, (b) (Subscript[S, z]), (c) parity or the total momentum, as well as (d) the total spin -- this is harder.

Hubbard[L_?IntegerQ,U_]:=Block[{lis={},N=4^L,v,v2},
lis=Table[{v,v}->U CountSinglets[v-1,L],{v,1,N}]; (* potential energy *)
For[v=0,v<N,v++, (* loop over initial state *)
For[n=1,n<L,n++, (* original site for hopping right *)
If[(v2=HopRight[v,n,2])≠-1,
lis=Append[lis,{v+1,v2+1}->- 1.0FermiFactorUpRight[v,n]]];
If[(v2=HopRight[v,n,1])≠-1,
lis=Append[lis,{v+1,v2+1}->- 1.0FermiFactorDnRight[v,n]]]];
For[n=0,n<L-1,n++, (* original site for hopping left *)
If[(v2=HopLeft[v,n,2])≠-1,
lis=Append[lis,{v+1,v2+1}->- 1.0FermiFactorUpLeft[v,n]]];
If[(v2=HopLeft[v,n,1])≠-1,
lis=Append[lis,{v+1,v2+1}->- 1.0FermiFactorDnLeft[v,n]]]];
];
(* can further optimize by using {v2,v} with left hop *)
SparseArray[Sort[lis],{N,N}]
]

HH=Hubbard[2,2.0];

HH

HH//MatrixForm

(HH-Transpose[HH])//Norm

I am not sure why this does not work

HH1=Developer`ToPackedArray[N[HH]]

Developer`PackedArrayQ[HH1]

H1=Developer`ToPackedArray[N[HH],Real]

Developer`PackedArrayQ[H1]

HH=Hubbard[5,2.0];

HH

(energ=Eigenvalues[HH])//Timing

energ//Sort

Eigenvalues[HH,15]//Timing

Eigenvalues[HH,-15]//Timing

Take[energ//Sort,15]

Take[energ//Sort,-15]

This is the way to get the ground state energy

Eigenvalues[HH,5,Method->{Arnoldi,Shift->-5.0000}]//Timing

HH=Hubbard[6,2.0];

Eigenvalues[HH,5,Method->{Arnoldi,Shift->-15.0000}]//Timing

Eigenvalues[HH,5,Method->{Arnoldi,Shift->-5.1000}]//Timing

Eigenvalues[HH,15,Method->{Arnoldi,Shift->-5.1000}]//Timing

Fixed-N, Fixed-Sz computation

Easiest way is to restrict initial states.  Make functions for counting spin-up and spin-down electrons

PP[135]

DigitCount[BitAnd[135,BitShiftLeft[1,2*5]-1],4,{2,3}]

CountUp[v_,n_]:=Apply[Plus,DigitCount[BitAnd[v,BitShiftLeft[1,2n]-1],4,{2,3}]]

CountDn[v_,n_]:=Apply[Plus,DigitCount[BitAnd[v,BitShiftLeft[1,2n]-1],4,{1,3}]]

For[v=0,v<31,v++,Print[PP[v],"  ",CountUp[v,5],"  ",CountDn[v,5]]]

Try to make program for diagonalizing Hubbard model with N electrons and z-component of spin S

Position[{1,2,3,4},3]

Hubbard[Nel_?IntegerQ, Sz_?IntegerQ, L_?IntegerQ,U_]:=Block[{lis={},N=4^L,n1,n2,v,v2,states,pos,pos2},
(* generate sorted list of states *)
For[v=0;len=0;states={},v<N,v++, (* loop over initial state *)
n1=CountDn[v,L];
n2=CountUp[v,L];
If[n1+n2==Nel && n2-n1==Sz,(* Ok state *)
(* state is not there yet *)
states=Append[states,v];
len++;
(* Print[PP[v]]; *)
]];
states=Sort[states];
For[v=0,v<N,v++, (* loop over initial state *)
n1=CountDn[v,L];
n2=CountUp[v,L];
If[n1+n2==Nel && n2-n1==Sz, (* Ok state *)
pos=Position[states,v][[1]][[1]];
lis=Append[lis,{pos,pos}->U CountSinglets[v,L]]; (* potential energy *)
For[n=1,n<L,n++, (* original site for hopping right *)
If[(v2=HopRight[v,n,2])≠-1,
pos2=Position[states,v2] [[1]][[1]];
lis=Append[lis,{pos,pos2}->- 1.0FermiFactorUpRight[v,n]]];
If[(v2=HopRight[v,n,1])≠-1,
pos2=Position[states,v2] [[1]][[1]];
lis=Append[lis,{pos,pos2}->- 1.0FermiFactorDnRight[v,n]]]];
For[n=0,n<L-1,n++, (* original site for hopping left *)
If[(v2=HopLeft[v,n,2])≠-1,
pos2=Position[states,v2] [[1]][[1]];
lis=Append[lis,{pos,pos2}->- 1.0FermiFactorUpLeft[v,n]]];
If[(v2=HopLeft[v,n,1])≠-1,
pos2=Position[states,v2] [[1]][[1]];
lis=Append[lis,{pos,pos2}->- 1.0FermiFactorDnLeft[v,n]]]];
]
];
(* can further optimize by using {pos2,pos} with left hop *)
(* can further optimize by first generating the sorted list of states and then running the loop -- faster search *)
SparseArray[Sort[lis],{len,len}]
]

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

HH

HH=Hubbard[1,1,8,0.0];HH//MatrixForms

Eigenvalues[HH]//Sort

-2Cos[ Pi/9]//N

Check that all is working right by comparing with non-interacting spectrum at U=0

HH=Hubbard[2,2,6,0.0]; Take[Eigenvalues[HH]//Sort,4]

-2Cos[ Pi/7]-2Cos[ Pi/7]//N

-2Cos[ Pi/7]-2Cos[ 2 Pi/7]//N

-2Cos[2 Pi/7]-2Cos[ 2 Pi/7]//N

-2Cos[ Pi/7]-2Cos[ 3 Pi/7]//N

HH=Hubbard[4,4,8,4.0];

8!/(4!)^2

70

HH

Eigenvalues[HH]//Sort

-2Cos[ Pi/9]-2Cos[ 2 Pi/9]-2Cos[ 3 Pi/9]-2Cos[ 4Pi/9]//N

HH=Hubbard[8,8,8,4.0];

HH

HH//MatrixForm

HH=Hubbard[6,4,5,4.0];

HH

HH//MatrixForm