Heisenberg model, complete diagonalization
Spin ring:  N=8, spin=1/2

©  Prof. Dr. Jürgen Schnack, Universität Bielefeld, 03. 12. 2019

Initialization

Set -2 J = 1; (AF)
Basis is one position longer to simulate NS+1=1.

NS=8;
spin=(1/2);
DimH=(((2*spin+1))^NS)

Basis=Table[{0,0,0,0,0,0,0,0,0},{n,1,DimH}]; (* an NS anpassen *)

Create basis (position N+1 is a repetition of N=1)

For[i=2,i<DimH+1,i++,
For[j=1,j<NS+1,j++,
Basis[[i,j]]=Basis[[i-1,j]]];
Basis[[i,1]]=Basis[[i,1]]+1;
j=1;
While[j<NS,
If[Basis[[i,j]]==(2*spin+1),Basis[[i,j]]=0;Basis[[i,j+1]]=Basis[[i,j+1]]+1;j++,j=NS]]
];
For[i=1,i<DimH+1,i++,
Basis[[i,NS+1]]=Basis[[i,1]]];

MatrixForm[Basis]

Create full Hamilton matrix

Diagonal = Ising model

HMatrix=Table[0.0,{i,1,DimH},{j,1,DimH}];
For[i=1,i<DimH+1,i++,
For[j=1,j<DimH+1,j++,
		HMatrix[[i,j]]=((∏_(m=1)^NS)(KroneckerDelta[Basis[[j,m]]-Basis[[i,m]]]))*
			(Sum[(spin-Basis[[j,l]])*(spin-Basis[[j,l+1]]),{l,1,NS}])
]];

Off-diagonal = XY model

For[i=1,i<DimH+1,i++,
For[j=1,j<DimH+1,j++,
For[k=1,k<NS+1,k++,
HilfsElement=Table[Basis[[j,n]],{n,1,NS+1}];
HilfsElement[[k]]=HilfsElement[[k]]-1;
HilfsElement[[k+1]]=HilfsElement[[k+1]]+1;
If[k==1,HilfsElement[[NS+1]]=HilfsElement[[1]]];
If[k==NS,HilfsElement[[1]]=HilfsElement[[NS+1]]];
HMatrix[[i,j]]=HMatrix[[i,j]]+0.5*
Product[KroneckerDelta[Basis[[i,m]]-HilfsElement[[m]]],{m,1,NS}]*
(Sqrt[(Basis[[j,k]])*(2*spin-Basis[[j,k]]+1)])*
(Sqrt[(2*spin-Basis[[j,k+1]])*(Basis[[j,k+1]]+1)])
]]];

For[i=1,i<DimH+1,i++,
For[j=1,j<DimH+1,j++,
For[k=1,k<NS+1,k++,
HilfsElement=Table[Basis[[j,n]],{n,1,NS+1}];
HilfsElement[[k]]=HilfsElement[[k]]+1;
HilfsElement[[k+1]]=HilfsElement[[k+1]]-1;
If[k==1,HilfsElement[[NS+1]]=HilfsElement[[1]]];
If[k==NS,HilfsElement[[1]]=HilfsElement[[NS+1]]];
HMatrix[[i,j]]=HMatrix[[i,j]]+0.5*
Product[KroneckerDelta[Basis[[i,m]]-HilfsElement[[m]]],{m,1,NS}]*
(Sqrt[(Basis[[j,k+1]])*(2*spin-Basis[[j,k+1]]+1)])*
(Sqrt[(2*spin-Basis[[j,k]])*(Basis[[j,k]]+1)])
]]];

If you are desperate, have a look at the full matrix!

(* MatrixForm[HMatrix] *)

Eigenvalues of H

EV=Sort[Chop[Eigenvalues[HMatrix]]]

Lanczos

T-Matrix

NL=10;
TMatrix=Table[0.0,{i,1,NL},{j,1,NL}];

Startvektor = 1. Vektor

r1=Table[RandomReal[]-0.5,{j,1,DimH}];
n1=r1.r1;
r1=r1/(Sqrt[n1]);
r1.r1

2. Vektor

TMatrix[[1,1]]=r1.HMatrix.r1;
r2=HMatrix.r1-r1*TMatrix[[1,1]];
n2=r2.r2;
r2=r2/(Sqrt[n2]);
r2.r2

r2.r1

TMatrix[[1,2]]=r1.HMatrix.r2;
TMatrix[[2,1]]=TMatrix[[1,2]];
TMatrix[[2,2]]=r2.HMatrix.r2;

k. Vektor

For[k=3,k<NL+1,k++,
rk=HMatrix.r2-r1*TMatrix[[k-2,k-1]]-r2*TMatrix[[k-1,k-1]];
rk normieren;
TMatrix[[k,k]]=;
r1=r2;
r2=rk;
Print[Sort[Chop[Eigenvalues[TMatrix]]]]];