(** INPUT **)

n=12; (* lattice size *)
J=-1.0; (* J > 0 - Antiferromagnetic, J < 0 - Ferromagnetic case *)
α=1.0; (* J anisotropy : α=1 -> Heisenberg limit, α=0-> Ising limit *)

(** INPUT END **)

(** BASIS CONSTRUCTION **)

nst=2^n; (* number of states *)
magnst=Table[{nu-n/2,Binomial[n,nu]},{nu,0,n}]; (* set of pairs {magnetization 'm', number of states for 'm'} *)
base=Table[{},Length[magnst]] ; (* base init - list of the same size as 'magnst' containing empty lists *)
Getnu[st_]:=(
nu=0;
For[k=0,k<BitLength[st],k++,If[BitGet[st,k]==1,nu++]];
nu); (* finds number of spins up 'nu' for state 'st' *)
For[st=0,st<nst,st++,AppendTo[base[[1+Getnu[st]]],st]]; (* ordering and splitting the base respectively to magnetization - magnetization conservation *)

(** BASIS CONSTRUCTION END **)

(** BLOCK HAMILTONIAN CONSTRUCTION **)

Needs["Combinatorica`"]; (* needed for BinarySearch procedure *)
BlockConstruct[m_]:=(
block=SparseArray[Table[0.0,{i,1,magnst[[m,2]]},{j,1,magnst[[m,2]]}]]; (* block init - matrix filled with 0 values *)
For[imst=1,imst<=magnst[[m,2]],imst++,( (* For each state indexed with 'imst' at given subspace 'm' ... *)
ist=base[[m,imst]];(* ... take its decimal index 'ist' (decimal representation) in binary base. *)
For[i=0,i<n,i++,( (* Loop over lattice sites 'i' from 0 to n-1 corresponding to (i^th) bit in binary representation of 'ist'. *)
j=Mod[i+1,n]; (* Get nearest neighbour (higher bit) 'j' of 'i' with periodic boundary conditions. *)
If[BitGet[ist,i]==BitGet[ist,j],block[[imst,imst]]+=0.25 J,( (* If spins (bits) at i and j are the same, add (J/4) to (Subscript[H, imst,imst]) matrix element ... *)
block[[imst,imst]]-=0.25 J; (* ... otherwise add -(J/4) to (Subscript[H, imst,imst]) matrix element, ... *)
jst=BitXor[ist,2^i+2^j]; (* ... reverse spins at 'i' and 'j' in state 'ist' to get state 'jst' ... *)
jmst=BinarySearch[base[[m]],jst]; (* ... and find its index 'jmst' in basis spanning subspace m. *)
block[[imst,jmst]]+=0.5 α J; (* Add ((Jα)/2) to (Subscript[H, imst,jmst]) matrix element. *)
)];
)];
)];
block); (* constructs a single block of the block Hamiltonian for a given subspace (magnetization), m = 1, 2, 3, ..., Length[magnst] *)
Hamiltonian={};
For[m=1,m≤ Length[magnst],m++,AppendTo[Hamiltonian, BlockConstruct[m]]] (* Loop over subspaces m to get full block Hamiltonian *)

(** BLOCK HAMILTONIAN CONSTRUCTION END **)

(** DIAGONALIZATION **)

DiagEneVec={}; (* eigensystem container - init *)
For[m=1,m≤ Length[magnst],m++,AppendTo[DiagEneVec,Eigensystem[Hamiltonian[[m,All]],Method->"Direct"]]] (* diagonalizes blocks separately *)
DiagHam=SparseArray[DiagonalMatrix[Level[DiagEneVec[[All,1]],{-1}]]]; (* combines diagonalized blocks into full diagonal Hamiltonian *)

(** DIAGONALIZATION END **)

(** GROUND ENERGY AND VECTORS **)

GSE=DiagEneVec[[1,1,1]]; (* ground state energy - init *)
GSV={}; (* ground state eigenvectors container -init *)
For[m=1,m≤ Length[magnst],m++, (* Within each subspace 'm' ...*)
For[mst=1,mst≤magnst[[m,2]],mst++, (* ... for each eigenstate 'mst' ... *)
EV=DiagEneVec[[m,1,mst]]; (* ... get its eigenvalue 'EV'. *)
If[EV<GSE,( (* If 'EV' is lower than assumed ground state energy 'GSE' ... *)
GSE=EV; (* ... then assume 'GSE' = 'EV', ... *)
GSV={{m,DiagEneVec[[m,2,mst]]}}; (* ... and assume new ground state eigenvector 'mst'. *)
),
If[EV==GSE,AppendTo[GSV,{m,DiagEneVec[[m,2,mst]]}]] (* Else if 'EV' = 'GSE' then append next ground state eigenvector 'mst' to 'GSV' *)
];
];
];

(** GROUND ENERGY AND VECTORS END **)

(** OPERATORS DEFINITION **)

Sz[i_,st_]:={BitGet[st,i]-0.5,st}; (* (Subscript[S, i]^z)(Ket[st]) : returns {±(1/2),(Ket[st])} depending on spin value at site i *)
Sp[i_,st_]:={1.0-BitGet[st,i],BitSet[st,i]};(* (Subscript[S, i]^(+))(Ket[st]) : returns {1,(Ket[st+(2^i)])} if spin at site i was down or returns {0,(Ket[st])} otherwise | (Ket[st+(2^i)])-(Ket[st])=(Ket[1 at site i and 0 everywhere else in bit representation]) *)
Sm[i_,st_]:={0.0+BitGet[st,i],BitClear[st,i]}; (* (Subscript[S, i]^(-))(Ket[st]) : returns {1,(Ket[st-(2^i)])} if spin at site i was up or returns {0,(Ket[st])} otherwise | (Ket[st])-(Ket[st-(2^i)])=(Ket[1 at site i and 0 everywhere else in bit representation]) *)

(* Sites numeration is i = 0, 1, ..., nst-1 *)

Sqp[q_,st_]:=Table[(s=Sp[r,st];{Exp[2 Pi I q r/n]s[[1]],s[[2]]}),{r,0,n-1}];
Sqm[q_,st_]:=Table[(s=Sm[r,st];{Exp[2 Pi I q r/n]s[[1]],s[[2]]}),{r,0,n-1}];

(** OPERATORS DEFINITION END **)

(** GREEN'S FUNCTION CALCULATIONS **)

G[k_,ω_,GS_]:=( (* 'k' & 'ω' -> tables of values for which Green's function will be calculated, 'GS' -> for example 'GSV[[1]]' *)
q=Mod[Round[n k / (2 Pi)],n];
m=GS[[1]]; (* magnetization of ground state -> see structure of 'GSV[[1]]' *)
mu=If[m==Length[magnst],m,m+1]; (* magnetization of a state obtained applying S+ to the ground state *)
mnst=magnst[[m,2]]; (* size (number of states) of the subspace with magnetization 'm' *)
munst=magnst[[mu,2]]; (* size (number of states) of the subspace with magnetization 'mu' *)

ESqp[q_,est_]:=( (* assumes eigenstate est is from the subspace of magnetization 'm', applies S+, result is a vector from the subspace 'mu' written in the binary basis of subspace 'mu' *)
SqpEst=Table[0.,munst];
For[mst=1,mst≤mnst,mst++,(
bst=base[[m,mst]];
SqpBst=Sqp[q,bst];
For[it=1,it≤n,it++,(
If[SqpBst[[it,1]]≠0.,(
must=BinarySearch[base[[mu]],SqpBst[[it,2]]];
SqpEst[[must]]+=est[[2,mst]]*SqpBst[[it,1]];
)];
)];
)];
SqpEst);

nw=Length[ω];
nq=Length[q];
result=Table[0.0,nw,nq];
denom=Table[ω[[iw]]+I δ-(DiagEneVec[[mu,1,vst]]-GSE),{iw,1,nw},{vst,1,munst}]; (* denominator does not depend on 'q' therefore it is calculated here for speed up *)
For[iq=1,iq≤nq,iq++,( (* loop over 'q' ... *)
SqpGS=ESqp[q[[iq]],GS];
For[vst=1,vst≤munst,vst++,( (* ... loop over eigen values and vectors ... *)
numer=Abs[Total[Conjugate[DiagEneVec[[mu,2,vst]]]*5SqpGS]]^2;
For[iw=1,iw≤nw,iw++,( (* ... loop over 'ω' *)
result[[iw,iq]]+=numer/denom[[iw,vst]];
)];
)];
)];
result);

(** GREEN'S FUNCTION CALCULATIONS END **)

(** DIAGONALIZATION RESULTS **)

Print["System size , Total number of states    ", n, " , ", nst];
Print[""];Print["Possible magnetizations    ", magnst[[All,1]]];
Print["Number of states for corresponding magnetization    ", magnst[[All,2]]];
Print[""];Print["Ground state energy , degeneracy    ", GSE, " , ", Length[GSV]];
Print[""];Print["Blocks of the Hamiltonian ordered from lowest to highest magnetization    "];
Table[ArrayPlot[Hamiltonian[[mag,All]]],{mag,1,Length[magnst]}]
Print[""];Print["Blocks of the Hamiltonian ordered from lowest to highest magnetization after diagonalization    "];
Table[ArrayPlot[DiagonalMatrix[Level[DiagEneVec[[mag,1]],{-1}]]],{mag,1,Length[magnst]}]
Print[""];Print["Full Hamiltonian after diagonalization of separate blocks    "];
Show[ArrayPlot[DiagHam],ImageSize->600]
(*Print[""];Print["Ground State    "];
Grid[{Table[Grid[{{Grid[{{"Magnetization: ",magnst[[GSV[[i,1]],1]]}}]},{Show[ArrayPlot[{GSV[[i,2]]}],ImageSize->50 n]}}],{i,1,Length[GSV]}]}]
 Grid[{{MatrixForm[GSV[[1,2]]], MatrixForm[base[[GSV[[1,1]]]]]}}]*)


(** DIAGONALIZATION RESULTS END **)

(** GREEN'S FUNCTION RESULTS **)

δ=0.0002; (* peak broadening *)
dk=1;
dw=0.1;
wl=-0.5;
wr=3.5;

ω=Table[iw,{iw,wl,wr,dw}];
k=Table[Pi(2 ik/n-1),{ik,0,n+1,dk}];

A=-1/Pi Im[G[k,ω,GSV[[1]]]];

ListDensityPlot[A,DataRange->{{-Pi- Pi/n,Pi + Pi/n},{wl-dw/2,wr-dw/2}},PlotRange->{0,1},InterpolationOrder->0,ClippingStyle->Automatic,Mesh->None,ColorFunction->"FuchsiaTones",ImageSize->800,AspectRatio->1/N[GoldenRatio]]

(** GREEN'S FUNCTION RESULTS END **)

(* Kiedy wyjaz? *)