Definions of the Rauzy sequence (to be run once for all)

ClearAll[R]

SetDirectory[NotebookDirectory[]];

R[-1]=0;
R[0]=1;
R[1]=1;

R[n_]:=R[n-1]+R[n-2]+R[n-3]

ClearAll[Rp,Rm]

Rp[-2]=0;
Rp[-1]=0;
Rp[0]=1;
Rp[n_]:=Rp[n-1]+Rp[n-2]+Rp[n-3]

Rm[0]=1;
Rm[-1]=1;
Rm[-2]=2;
Rm[n_]:=Rm[n-3]-Rm[n-2]-Rm[n-1]

Hamiltonian with PBC

order=11;
q=R[order-2];
r=R[order-1];
s=R[order];
ntot=R[order+1]

Spectrum in zero field

c=0;
rawham=Reap[
Do[c=c+1;
Sow[{{i,i+q},1}];
Sow[{{i+q,i},1}],
{i,1,ntot-q}
]

Do[c=c+1;
Sow[{{i,i+r},1}];
Sow[{{i+r,i},1}],
{i,1,ntot-r}
]
Do[c=c+1;
Sow[{{i,i+s},1}];
Sow[{{i+s,i},1}],
{i,1,ntot-s}
]

][[2,1]];

tmp=Split[Sort[rawham],#1[[1]]==#2[[1]]&];
rawham=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];

ham=SparseArray[rawham];
Length[ham]

res=Sort[Eigenvalues[N[-ham]]];
Export["spec_"<>ToString[order]<>"_pbc_0_0_1.dat",res]

res[[1]]

nb=50;
Print["Symmetric spectrum : ",Union[Chop[res+Reverse[res]]]=={0}]
Histogram[Chop[res],{Delete[Table[4.2i/(nb),{i,-nb,nb}],{nb+1}]}]

Variational ansatz 1

Finite-size variational state (vecana=vec)

tab=Table[Tr[ham[[i]]],{i,1,ntot}];
Count[tab,3]==2q
Count[tab,4]==ntot-4q
Count[tab,5]==2q
vec=Table[Switch[tab[[i]],3,α,4,β,5,γ],{i,1,ntot}];
Normvec=Simplify[Norm[vec],{α∈Reals,β∈Reals,γ∈Reals}];
vec=vec/Normvec;

vecana=Join[
Table[α,{i,1,q}],
Table[β,{i,q+1,r}],
Table[γ,{i,r+1,ntot-s}],
Table[β,{i,ntot-s+1,s}],
Table[γ,{i,s+1,ntot-r}],
Table[β,{i,ntot-r+1,ntot-q}],
Table[α,{i,ntot-q+1,ntot}]];
test=vecana;
Normvecana=(Sqrt[2q (α^2)+(ntot-4q)(β^2)+2q (γ^2)]);
vecana=vecana/Normvecana;
Union[Simplify[vec-vecana]]=={0}

hamtest=Join[Table[β+2γ,{i,1,r-q}],
Table[3γ,{i,r-q+1,q}],
Table[α+2 β+γ,{i,q+1,s-r}],
Table[α+β+2 γ,{i,s-r+1,r}],
Table[3 α+β+γ,{i,r+1,2q}],
Table[2 α+2 β+γ,{i,2q+1,s-q}],
Table[2 α+ β+2γ,{i,s-q+1,ntot-2r}],
Table[2 α+2 β+γ,{i,ntot-2r+1,ntot-s}],

Table[2 β+2 γ,{i,ntot-s+1,s}],

Table[2 α+ 2β+γ,{i,s+1,2r}],
Table[2 α+ β+2γ,{i,2r+1,ntot-(s-q)}],
Table[2 α+2 β+γ,{i,ntot-(s-q)+1,ntot-2q}],
Table[3 α+β+γ,{i,ntot-2q+1,ntot-r}],
Table[α+β+2 γ,{i,ntot-r+1,ntot-(s-r)}],
Table[α+2β+ γ,{i,ntot-(s-r)+1,ntot-q}],
Table[3γ,{i,ntot-q+1,ntot-(r-q)}],
Table[β+2γ,{i,ntot-(r-q)+1,ntot}]
]==ham.test

Tr[ham.ham]/ntot

Finite-size variational energy Evar = vec.ham.vac

Evar=(2(r-q)α(β+2γ)+2 (2q-r)3α γ+2β(s-r-q)(α+2 β+γ)+2β(2r-s)(α+β+2 γ)+2γ(2q-r)(3 α+β+γ)+2γ(s-3q+2r-s)(2 α+2 β+γ)+2γ(ntot-2r-s+q)(2 α+ β+2γ)+β(2s-ntot)(2β+2γ))/(2q (α^2)+(ntot-4q)(β^2)+2q (γ^2)) /.ntot->q+r+s;

Evarsimp=(((6q-2r)((γ^2)-(β^2))+4(s-r) (β^2)+4(r-q)β(α+2γ)+4α γ(4q-r))/(2 q ((α^2)+(γ^2))+(-3 q+r+s) (β^2)));
Simplify[Evar-Evarsimp]

Simplify[Evar-vec.ham.vec]

NMinimize[-Evarsimp,{α,β,γ}]

Variational energy in the thermodynamical limit

(*Divide by q for thermodynamical limit + r/q->θ/.s/q->(θ^2) where (θ^3)=(θ^2)+θ+1 *)
sol=Solve[(θ^3)-(θ^2)-θ-1==0,θ]
Evarsimpthermo=(( (β^2) (4(θ^2)-2θ-6)+(γ^2) (6-2 θ)+4 α γ (4-θ)+4 β (α+2 γ) (θ-1))/(2 ((α^2)+(γ^2))+(β^2) ((θ^2)+θ-3)))/.sol[[1]];
sol2=NMinimize[-Evarsimpthermo/.sol[[1]],{α,β,γ}]

Overlap between the exact finite-size ground state (vecnum) and the finite-size variational ansatz (vecanamin) versus 1/order

vecanamin=vecana/.NMinimize[-Evar,{α,β,γ}][[2]];
vecnum=Eigensystem[N[-ham],10,Method->{"Arnoldi",MaxIterations->∞}][[2,1]];
Abs[vecnum.vecanamin]

listfid={{1/8,0.9997792457437525`},{1/9,0.9997781981592854`},{1/10,0.9997770826046033`},{1/11,0.9997760357687865`},{1/12,0.9997754947458526`},{1/13,0.9997754637524254`},{1/14,0.9997753404590759`},{1/15,0.9997751144826738`  },{1/16, 0.9997751636605063`  }} ;
ListPlot[listfid]

ListPlot[{vecnum,vecanamin},Joined->True]

Exact finite-size ground state energy versus 1/order

(*order=6;
q=R[order-2];
r=R[order-1];
s=R[order];
ntot=R[order+1]

c=0;
rawham=Reap[
Do[c=c+1;
Sow[{{i,i+q},1}];
Sow[{{i+q,i},1}],
{i,1,ntot-q}
]

Do[c=c+1;
Sow[{{i,i+r},1}];
Sow[{{i+r,i},1}],
{i,1,ntot-r}
]
Do[c=c+1;
Sow[{{i,i+s},1}];
Sow[{{i+s,i},1}],
{i,1,ntot-s}
]

(* Trick to shift the ground-state energy as the absolute minimum energy *)

Do[
Sow[{{i,i},10}],
{i,1,ntot}
]

][[2,1]];

tmp=Split[Sort[rawham],#1[[1]]==#2[[1]]&];
rawham=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];

ham=SparseArray[rawham];
res=Eigenvalues[N[-ham],100,Method -> {"Arnoldi",MaxIterations->50000}]+10;
res[[1]]*)

gsePBC={{1/12,-4.115003008951662`},{1/13,-4.11500907},{1/15,-4.115008095924509`},{1/16,-4.115008539510098` },{1/17,-4.115008493452219`},{1/18,-4.115008443234},{1/19,-4.115008473560783` },{1/20,-4.115008319349304`}}

ListPlot[gsePBC,PlotRange->{-4.115008,-4.1150091},Frame->True]

Variational ansatz 2 (second shell)

ham2=Table[ham[[i,j]],{i,1,ntot},{j,1,ntot}];

bla=Reap[Do[a=Tr[ham2[[i]]];
listpos=Position[ham2[[i]],1];
tmp={};
Do[AppendTo[tmp,Tr[ham2[[listpos[[j,1]]]]]],{j,1,Length[listpos]}];
Sow[{a,tmp}];
,{i,1,ntot}]][[2,1]];

Length[Union[bla]]
Length[Union[Sort[bla]]]

Clear[list,res,a]
list=Union[bla]
res=bla;
Do[res=res/.list[[i]]->a[i],{i,1,Length[list]}]

vecana2=res/Norm[res];
nrj=-vecana2.ham.vecana2;
NMinimize[nrj,Table[a[i],{i,1,Length[list]}]]

v1=-4.112973010832521`
v2=-4.114927329157021`
v0=-4.1150084

(v1-v0)/v0
(v2-v0)/v0

Overlap between the exact finite-size ground state (vecnum) and the finite-size variational ansatz (vecanamin) versus 1/order

vecanamin2=vecana2/.NMinimize[nrj,Table[a[i],{i,1,Length[list]}]][[2]];
vecnum=Eigensystem[N[-ham],10,Method->{"Arnoldi",MaxIterations->∞}][[2,1]];
Abs[vecnum.vecanamin2]

ListPlot[{vecnum,vecanamin,vecanamin2},PlotStyle->{Red,Green,Blue},Joined->True]

Variational ansatz 3 (third shell)

Clear[a,b,c,tmp1,tmp2]

bla=Reap[Do[
a=Tr[ham2[[i]]];
listpos1=Position[ham2[[i]],1];
tmp1={};
tmp2={};
Do[b=Tr[ham2[[listpos1[[j,1]]]]];
listpos2=Position[ham2[[listpos1[[j,1]]]],1];

AppendTo[tmp1,b];

Do[c=Tr[ham2[[listpos2[[k,1]]]]];
AppendTo[tmp2,c];
,{k,1,Length[listpos2]}];

,{j,1,Length[listpos1]}];

Sow[{a,tmp1,tmp2}];
,{i,1,ntot}]][[2,1]];

Length[Union[bla]]

ClearAll[list,res,a]
list=Union[bla]
Length[list]
res=bla;
Do[res=res/.list[[i]]->a[i],{i,1,Length[list]}];

Do[
Do[
If[Sort[list[[i,3]]]==Sort[list[[j,3]]],{Print[i," ",j],a[j]=a[i]}],{j,i+1,Length[list]}],
{i,1,Length[list]}]

vecana3=res/Norm[res];
nrj=-vecana3.ham.vecana3;
NMinimize[nrj,Table[a[i],{i,1,Length[list]}]]

Overlap between the exact finite-size ground state (vecnum) and the finite-size variational ansatz (vecanamin) versus 1/order

vecanamin3=vecana3/.NMinimize[nrj,Table[a[i],{i,1,Length[list]}]][[2]];
vecnum=Eigensystem[N[-ham],50,Method->{"Arnoldi",MaxIterations->∞}][[2,1]];
Abs[vecnum.vecanamin3]

ListPlot[Abs[{vecnum,vecanamin,vecanamin2,vecanamin3}],PlotStyle->{Red,Green,Blue,Magenta},Joined->True]

Clear[a,res,list];
bla=Table[MatrixPower[ham,10][[i,i]],{i,1,ntot}];
list=Union[bla];
Length[list]
res=bla;
Do[res=res/.list[[i]]->a[i],{i,1,Length[list]}]
vectest=res/Norm[res];

nrj=-vectest.ham.vectest;

NMinimize[nrj,Table[a[i],{i,1,Length[list]}]][[1]]

bla2=Table[MatrixPower[ham,2][[i,i]],{i,1,ntot}]
bla4=Table[MatrixPower[ham,4][[i,i]],{i,1,ntot}]
bla6=Table[MatrixPower[ham,6][[i,i]],{i,1,ntot}]
bla6=Table[MatrixPower[ham,8][[i,i]],{i,1,ntot}]

Position[bla2,4]
Position[bla4,27]
Position[bla6,315]

!

Position[bla6,315]

bla4[[68]]
bla6[[13]]
bla6[[24]]

Hamiltonian with OBC

order=17;
q=R[order-2];
r=R[order-1];
s=R[order];
ntot=R[order+1]

Construction

ClearAll[q,r,s,ntot,Aperp,A1,A2,gene,coord,coordiso]

q=Rp[order-2];
r=Rp[order-1];
s=Rp[order];
ntot=Rp[order+1]

Aperp={Rp[order],Rp[order-1],Rp[order-2]};
gene={Rm[order],Rm[order+1]+Rm[order+2],Rm[order+1]};
A1=    {Rm[order+1],Rm[order+2]+Rm[order+3],Rm[order+2]};
A2=    {Rm[order+2],Rm[order+3]+Rm[order+4],Rm[order+3]};

Det[{A1,A2,gene}]==1;

L=Transpose[{A1,A2,Aperp}];
Linv=Inverse[L];

a={1,0,(1/(Sqrt[3]))};
b={-(1/2),((Sqrt[3])/2),(1/(Sqrt[3]))};
c={-(1/2),-((Sqrt[3])/2),(1/(Sqrt[3]))};

M=Transpose[{a,b,c}];

coord[j_]:=Drop[FractionalPart[j*Linv.gene],-1];

coordiso[j_]:=Drop[M.L.FractionalPart[j*Linv.gene],-1];

A1iso=Drop[M.A1,-1];
A2iso=Drop[M.A2,-1];
geneiso=Drop[M.gene,-1];

k=1.;
Norm[coord[k]-coord[k+q]]
Norm[coord[k]-coord[k+r]]
Norm[coord[k]-coord[k+s]]

Norm[coordiso[k]-coordiso[k+q]]
Norm[coordiso[k]-coordiso[k+r]]
Norm[coordiso[k]-coordiso[k+s]]

Canonical basis (Subscript[A, 1])={1,0}, (Subscript[A, 2])={0,1}

ϵ=0.5;
graph=Reap[
Do[{
distance=Norm[coord[j]-coord[j+q]];
If[distance<ϵ,Sow[Line[{coord[j],coord[j+q]}]]];
},
{j,1,ntot-q}]

Do[{
distance=Norm[coord[j]-coord[j+r]];
If[distance<ϵ,Sow[Line[{coord[j],coord[j+r]}]]];
},
{j,1,ntot-r}]

Do[{
distance=Norm[coord[j]-coord[j+s]];
If[distance<ϵ,Sow[Line[{coord[j],coord[j+s]}]]];
},
{j,1,ntot-s}]

][[2,1]];

Graphics[{graph,Line[{{0,0},{0,-1},{-1,-1},{-1,0},{0,0}}]}]

Isometric basis

ϵ=1.1;
shift=1;
graphiso=Reap[
Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+q-shift]];
If[distance<ϵ,Sow[Line[{coordiso[j-shift],coordiso[j+q-shift]}]]];
},
{j,1,ntot-q}]

Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+r-shift]];
If[distance<ϵ,Sow[Line[{coordiso[j-shift],coordiso[j+r-shift]}]]];
},
{j,1,ntot-r}]

 Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+s-shift]];
If[distance<ϵ,Sow[Line[{coordiso[j-shift],coordiso[j+s-shift]}]]];
},
{j,1,ntot-s}]

][[2,1]];

P1=Graphics[{Red,Line[{{0,0},{0,1}}],Line[{{0,0},{1,0}}]}];
P2=Graphics[{graphiso,Line[{{0,0},-A1iso,-A1iso-A2iso,-A2iso,{0,0}}]}];
Show[P1,P2,P1]

Spectrum in zero field

ϵ=1.2;
c=0;
shift=1;
f=0;
rawham=Reap[
Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+q-shift]];
If[distance<ϵ,Sow[{{j,j+q},Exp[-I(2*Pi*f/Sqrt[3]*(coordiso[j+q-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+q-shift][[2]]))]}];c=c+1;Sow[{{j+q,j},Exp[I(2*Pi*f/Sqrt[3]*(coordiso[j+q-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+q-shift][[2]]))]}]];
},
{j,1,ntot-q}]

Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+r-shift]];
If[distance<ϵ,Sow[{{j,j+r},Exp[-I(2*Pi*f/Sqrt[3]*(coordiso[j+r-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+r-shift][[2]]))]}];c=c+1;Sow[{{j+r,j},Exp[I(2*Pi*f/Sqrt[3]*(coordiso[j+r-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+r-shift][[2]]))]}]];
},
{j,1,ntot-r}]

Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+s-shift]];
If[distance<ϵ,Sow[{{j,j+s},Exp[-I(2*Pi*f/Sqrt[3]*(coordiso[j+s-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+s-shift][[2]]))]}];c=c+1;Sow[{{j+s,j},Exp[I(2*Pi*f/Sqrt[3]*(coordiso[j+s-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+s-shift][[2]]))]}]];
},
{j,1,ntot-s}]

][[2,1]];

tmp=Split[Sort[rawham],#1[[1]]==#2[[1]]&];
rawham=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];
ham=SparseArray[rawham];

Length[ham]
res=Sort[Eigenvalues[N[ham]]];
ListPlot[res]
Print[c]

nb=51;
Print["Symmetric spectrum : ",Union[Chop[res+Reverse[res]]]=={0}]
Histogram[Chop[res],{Delete[Table[4.2i/(nb),{i,-nb,nb}],{nb+1}]}]

Tr[ham[[1]]]

Filling the matrix in a field

ϵ=1.01;
shift=1;

ClearAll[rawham,f];

rawham=Reap[
Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+q-shift]];
If[distance<ϵ,Sow[{{j,j+q},Exp[-I(2*π*f/(Sqrt[3])*(coordiso[j+q-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+q-shift][[2]]))]}];Sow[{{j+q,j},Exp[I(2*π*f/(Sqrt[3])*(coordiso[j+q-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+q-shift][[2]]))]}]];
},
{j,1,ntot-q}]

Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+r-shift]];
If[distance<ϵ,Sow[{{j,j+r},Exp[-I(2*π*f/(Sqrt[3])*(coordiso[j+r-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+r-shift][[2]]))]}];Sow[{{j+r,j},Exp[I(2*π*f/(Sqrt[3])*(coordiso[j+r-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+r-shift][[2]]))]}]];
},
{j,1,ntot-r}]

Do[{
distance=Norm[coordiso[j-shift]-coordiso[j+s-shift]];
If[distance<ϵ,Sow[{{j,j+s},Exp[-I(2*Pi*f/Sqrt[3]*(coordiso[j+s-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+s-shift][[2]]))]}];c=c+1;Sow[{{j+s,j},Exp[I(2*Pi*f/Sqrt[3]*(coordiso[j+s-shift][[1]]+coordiso[j-shift][[1]])*(coordiso[j-shift][[2]]-coordiso[j+s-shift][[2]]))]}]];
},
{j,1,ntot-s}]
][[2,1]];

tmp=Split[Sort[rawham],#1[[1]]==#2[[1]]&];
rawham=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];
ham=SparseArray[rawham];

Diagonalization in a field

nbpoints=256;
bc=obc;
fmin=0;
fmax=0.01;

data=Monitor[Reap[
For[i=1,i≤nbpoints+1,i++,
rulef={f->fmin+fmax(i-1.)/nbpoints};
ham=SparseArray[rawham/.rulef];
Sow[Flatten[{f,Sort[Eigenvalues[ham]]}]/.rulef]
(*
res=Eigensystem[ham];
For[j=1,j≤ntot,j++,
If[Tr[((Abs[res[[2,j]]])^4)]<ϵ,Sow[Flatten[{f,Chop[res[[1,j]]]}]/.rulef]];
];*)

]][[2,1]],i];

Dimensions[data]

Export["spec_"<>ToString[order]<>"_"<>ToString[bc]<>"_"<>ToString[fmin]<>"_"<>ToString[fmax]<>"_"<>ToString[nbpoints+1]<>".dat",data];

Landau level analysis

nbpoints=256;
bc=obc;
fmin=0;
fmax=0.01;

order=17;

bla=Import["spec_"<>ToString[order]<>"_"<>ToString[bc]<>"_"<>ToString[fmin]<>"_"<>ToString[fmax]<>"_"<>ToString[nbpoints+1]<>".dat","Table"];

tbp=Table[ListPlot[Transpose[{Transpose[bla][[1]],Transpose[bla][[i]]}],PlotStyle->{Blue,Thickness[0.0003]},PlotRange->{{0,0.01},{-4.12,-3.85}},Joined->True],{i,2,800}];

coeff=1.95735*2π/(Sqrt[3]) (*fit with the "exact" effective mass computed below *) 
papillon=Show[tbp,Plot[Table[{-4.1150+n coeff x},{n,1,37,2}],{x,0,0.01},PlotStyle->Red],ImageSize->800]

Clear[k]
(*fit with the effective mass *) 
kmin=10;
kmax=20;
tt=Take[bla,{kmin,kmax}];
list=Transpose[{Transpose[tt][[1]],Transpose[tt][[2]]}];
fit=Fit[list,{1,x},x]
coeff=Coefficient[fit,x]
Correlation[list]

Show[tbp,Plot[Table[{fit+(n-1) coeff x},{n,1,37,2}],{x,0,0.01},PlotStyle->Red],ImageSize->800]

Ground-state energy at f=0

nbpoints=512;
bc=obc;
fmin=0;
fmax=1;

listgse={};

Monitor[Do[
data=Import["gse_"<>ToString[order]<>"_"<>ToString[bc]<>"_"<>ToString[fmin]<>"_"<>ToString[fmax]<>"_"<>ToString[nbpoints+1]<>".dat","Table"];

ntot=Rp[order+1];
AppendTo[listgse,{1/order,data[[1,2]]}],
{order,10,18}]
,order]
listgse

ListPlot[listgse]

Ground-state energy vs flux

nbpoints=512;
bc=obc;
fmin=0;
fmax=1;

order=18;

data=Import["gse_"<>ToString[order]<>"_"<>ToString[bc]<>"_"<>ToString[fmin]<>"_"<>ToString[fmax]<>"_"<>ToString[nbpoints+1]<>".dat","Table"];tbp18=Table[ListPlot[Transpose[{Transpose[data][[1]],Transpose[data][[i]]}],PlotStyle->{Black,PointSize[0.0002]},PlotRange->{-4,-2.6},Joined->False],{i,2,2}];

Show[tbp18,ImageSize->800]

Effective mass tensor for the lowest mini-band (1/m=1.95735(1))

order=16;
q=R[order-2];
r=R[order-1];
s=R[order];
ntot=R[order+1]

Construction

ClearAll[q,r,s,ntot,Aperp,A1,A2,gene,coord,coordiso]

q=Rp[order-2];
r=Rp[order-1];
s=Rp[order];
ntot=Rp[order+1]

Aperp={Rp[order],Rp[order-1],Rp[order-2]};
gene={Rm[order],Rm[order+1]+Rm[order+2],Rm[order+1]};
A1=    {Rm[order+1],Rm[order+2]+Rm[order+3],Rm[order+2]};
A2=    {Rm[order+2],Rm[order+3]+Rm[order+4],Rm[order+3]};

Det[{A1,A2,gene}]==1;

L=Transpose[{A1,A2,Aperp}];
Linv=Inverse[L];

a={1,0,(1/(Sqrt[3]))};
b={-(1/2),((Sqrt[3])/2),(1/(Sqrt[3]))};
c={-(1/2),-((Sqrt[3])/2),(1/(Sqrt[3]))};

M=Transpose[{a,b,c}];

coord[j_]:=Drop[FractionalPart[j*Linv.gene],-1];

coordiso[j_]:=Drop[M.L.FractionalPart[j*Linv.gene],-1];

A1iso=Drop[M.A1,-1];
A2iso=Drop[M.A2,-1];
geneiso=Drop[M.gene,-1];

k=1.;
Norm[coord[k]-coord[k+q]]
Norm[coord[k]-coord[k+r]]
Norm[coord[k]-coord[k+s]]

Norm[coordiso[k]-coordiso[k+q]]
Norm[coordiso[k]-coordiso[k+r]]
Norm[coordiso[k]-coordiso[k+s]]

Scan kx, ky, kx=ky

Norm[coordiso[1+q]-coordiso[1]]
Norm[coordiso[1+r]-coordiso[1]]
Norm[coordiso[1+s]-coordiso[1]]
listx={};
nbpoints=10;
Monitor[Do[
kx=-Pi/(50Norm[A1])+(ii-1)/nbpoints*2 Pi/(50Norm[A1]);
k={kx,0};
aa=Exp[I k.(coordiso[1+q]-coordiso[1])];
bb=Exp[I k.(coordiso[1+r]-coordiso[1])];
cc=Exp[I k.(coordiso[1+s]-coordiso[1])];

rawham=Reap[
Do[Sow[{{j,j+q},aa}];Sow[{{j+q,j},Conjugate[aa]}],{j,1,ntot-q}];
Do[Sow[{{j,j+r},bb}];Sow[{{j+r,j},Conjugate[bb]}],{j,1,ntot-r}];
Do[Sow[{{j,j+s},cc}];Sow[{{j+s,j},Conjugate[cc]}],{j,1,ntot-s}];
][[2,1]];

tmp=Split[Sort[rawham],#1[[1]]==#2[[1]]&];
rawham=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];
ham=SparseArray[rawham];
res2=Sort[Eigenvalues[N[-ham],3,Method->{"Arnoldi",MaxIterations->20000}]];
AppendTo[listx,{kx,Min[res2]}];
,{ii,1,nbpoints}],ii]

listy={};
nbpoints=10;
Monitor[Do[
ky=-Pi/(50Norm[A2])+(ii-1)/nbpoints*2 Pi/(50Norm[A2]);
k={0,ky};
aa=Exp[I k.(coordiso[1+q]-coordiso[1])];
bb=Exp[I k.(coordiso[1+r]-coordiso[1])];
cc=Exp[I k.(coordiso[1+s]-coordiso[1])];

rawham=Reap[
Do[Sow[{{j,j+q},aa}];Sow[{{j+q,j},Conjugate[aa]}],{j,1,ntot-q}];
Do[Sow[{{j,j+r},bb}];Sow[{{j+r,j},Conjugate[bb]}],{j,1,ntot-r}];
Do[Sow[{{j,j+s},cc}];Sow[{{j+s,j},Conjugate[cc]}],{j,1,ntot-s}];
][[2,1]];

tmp=Split[Sort[rawham],#1[[1]]==#2[[1]]&];
rawham=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];
ham=SparseArray[rawham];
res2=Sort[Eigenvalues[N[-ham],3,Method->{"Arnoldi",MaxIterations->20000}]];
AppendTo[listy,{ky,Min[res2]}];
,{ii,1,nbpoints}],ii]

listxy={};
nbpoints=10;
Monitor[Do[
kxy=-Pi/(50Norm[A1])+(ii-1)/nbpoints*2 Pi/(50Norm[A1]);
k={kxy,kxy};
aa=Exp[I k.(coordiso[1+q]-coordiso[1])];
bb=Exp[I k.(coordiso[1+r]-coordiso[1])];
cc=Exp[I k.(coordiso[1+s]-coordiso[1])];

rawham=Reap[
Do[Sow[{{j,j+q},aa}];Sow[{{j+q,j},Conjugate[aa]}],{j,1,ntot-q}];
Do[Sow[{{j,j+r},bb}];Sow[{{j+r,j},Conjugate[bb]}],{j,1,ntot-r}];
Do[Sow[{{j,j+s},cc}];Sow[{{j+s,j},Conjugate[cc]}],{j,1,ntot-s}];
][[2,1]];

tmp=Split[Sort[rawham],#1[[1]]==#2[[1]]&];
rawham=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];
ham=SparseArray[rawham];
res2=Sort[Eigenvalues[N[-ham],3,Method->{"Arnoldi",MaxIterations->20000}]];
AppendTo[listxy,{kxy,Min[res2]}];
,{ii,1,nbpoints}],ii]

ListPlot[{listx,listxy,listy},PlotStyle->{Red,Green,Blue},PlotRange->{-4.12,-4.1}]

cx=Coefficient[Fit[listx,{1,x,x^2},x],x^2]
cy=Coefficient[Fit[listy,{1,x,x^2},x],x^2]
cxy=(Coefficient[Fit[listxy,{1,x,x^2},x],x^2]-cx-cy)/2
mat={{cx,cxy},{cxy,cy}};
Eigenvalues[mat]
2*(Sqrt[Det[mat]])

ListPlot[{{1/1705,1.9573424057149076` },{1/3706,1.9573514460669672` },{1/5768,1.9573507590182315` },{1/10609,1.9573471803025269` },{1/19513,1.9573486502404953` }}]

Orbital Susceptibility OBC

(* Ω[T_,i_,μ_]:=-(T/ntot) Sum[Log[1+Exp[-(1/T)(data[[i,j]]-μ)]],{j,2,ntot+1}]
χ[T_,μ_]:=Module[{tmp,fit},
tmp=Table[{(i-1)fmax/nbpoints,Ω[T,i,μ]},{i,1,imax-imin+1}];
fit=Fit[tmp,{1,x,(x^2)},x];
-2Coefficient[fit,(x^2)]]

(bpoints=16;
bc=obc;
fmin=0;
fmax=0.01;

temp=0.0125;
imin=4;
imax=8;           
δμ=0.01;        

Monitor[Do[

data=Take[Import["spec_"<>ToString[order]<>"_obc_"<>ToString[bc]<>"_"<>ToString[fmin]<>"_"<>ToString[fmax]<>"_"<>ToString[nbpoints+1]<>".dat","Table"],{imin,imax}];

ntot=Rp[order+1];

res=Monitor[Table[{μ,χ[temp,μ]},{μ,-5,+5,δμ}],μ];

Export["chi_orb_"<>ToString[order]<>"_"<>ToString[temp]<>".dat",res],

{order,10,18}],order]*)

colors={Red,Green,Blue,Magenta,Brown,Cyan,Black};
listorder={10,11,12,13,14,15,16,17,18};
listtemp={0.0125,0.025,0.05,0.1,0.2};
colors={ColorData["DeepSeaColors"][1-i*0.2]};

Do[
Do[
order=listorder[[i]];
temp=listtemp[[j]];
data=Import["chi_orb_"<>ToString[order]<>"_obc_"<>ToString[temp]<>".dat","Table"];
P[order,temp]=ListPlot[data,Joined->True,PlotStyle->ColorData["Rainbow"][(i-1)/(Length[listorder]-1)],PlotRange->All,ImageSize->700];
P2[order,temp]=ListPlot[data,Joined->True,PlotStyle->ColorData["Rainbow"][(j-1)/(Length[listtemp]-1)],PlotRange->All,ImageSize->700],

{j,1,Length[listtemp]}],
{i,1,Length[listorder]}]

temp=0.05;
tab={16,17,18}
bla=Show[Table[P[order,temp],{order,Flatten[tab]}],ImageSize->800,PlotRange->All]
(*Export["chi_0.2.pdf",bla]*)

bla=Show[Table[P2[18,listtemp[[i]]],{i,3,5}],ImageSize->800,PlotRange->{-6,All}]

Spin Susceptibility OBC

ClearAll[data,temp,μ,ntot];

χS[T_,μ_]:=Sum[(Sech[(data[[i]]-μ)/(2T)])^2,{i,2,ntot+1}]/(2*Sqrt[3]*T*ntot); 

nbpoints=16;
bc=obc;
fmin=0;
fmax=0.01;

temp=0.001;
δμ=0.01;        

Monitor[Do[

data=Import["spec_"<>ToString[order]<>"_"<>ToString[bc]<>"_"<>ToString[fmin]<>"_"<>ToString[fmax]<>"_"<>ToString[nbpoints+1]<>".dat","Table"][[1]];

ntot=Rp[order+1];

resfinal=Table[{μ,χS[temp,μ]},{μ,-5,+5,δμ}];

Export["chi_spin_"<>ToString[order]<>"_obc_"<>ToString[temp]<>".dat",resfinal],

{order,10,18}],{order,μ}]

listorder={10,11,12,13,14,15,16,17,18};
listtemp={0.001,0.0125,0.025,0.05,0.1,0.2};

Do[
Do[
order=listorder[[i]];
temp=listtemp[[j]];
data=Import["chi_spin_"<>ToString[order]<>"_obc_"<>ToString[temp]<>".dat","Table"];

Pobc[order,temp]=ListPlot[data,Joined->True,PlotStyle->ColorData["Rainbow"][(i-1)/(Length[listorder]-1)],PlotRange->All,ImageSize->700];

Pobc2[order,temp]=ListPlot[data,Joined->True,PlotStyle->ColorData["Rainbow"][(j-1)/(Length[listtemp]-1)],PlotRange->All,ImageSize->700],

{j,1,Length[listtemp]}],
{i,1,Length[listorder]}]

temp=0.0125;
tab={17,18}
bla=Show[Table[P[order,temp],{order,Flatten[tab]}],ImageSize->800,PlotRange->All]

bla=Show[Table[P2[18,listtemp[[j]]],{j,1,5}],ImageSize->800,PlotRange->{{-4.13,-4.1},{0.03,0.05}},Frame->True]

m=1/2.1;
T=0.025;blo=Plot[{m/(2 Pi),NIntegrate[1/(4T*2*Pi)*k*(Sech[k^2/(4 m T)-(mu+4.114458244825543`)/(2 T)])^2,{k,0,Infinity}]},{mu,-4.5,-3.8},PlotRange->{All,All},PlotStyle->Magenta];

NIntegrate[1/(4T*2*Pi)*k*(Sech[k^2/(4 m T)-(-4.11+4.1150)/(2 T)])^2,{k,0,Infinity}]

Spin Susceptibility PBC

ClearAll[data,T,temp,μ,ntot];

χS[T_,μ_]:=Sum[(Sech[(data[[i]]-μ)/(2T)])^2,{i,1,ntot}]/(2*Sqrt[3]*T*ntot); 

nbpoints=1;
bc=pbc;
fmin=0;
fmax=0;

temp=0.001;
δμ=0.01;        

Monitor[Do[

data=Flatten[Import["spec_"<>ToString[order]<>"_"<>ToString[bc]<>"_"<>ToString[fmin]<>"_"<>ToString[fmax]<>"_"<>ToString[nbpoints]<>".dat","Table"]];

ntot=Rp[order+1];

resfinal=Table[{μ,χS[temp,μ]},{μ,-5,+5,δμ}];

Export["chi_spin_"<>ToString[order]<>"_pbc_"<>ToString[temp]<>".dat",resfinal],

{order,10,18}],{order,μ}]

listorder={10,11,12,13,14,15,16,17,18};
listtemp={0.0125,0.025,0.05,0.1,0.2};

Do[
Do[
order=listorder[[i]];
temp=listtemp[[j]];
data=Import["chi_spin_"<>ToString[order]<>"_pbc_"<>ToString[temp]<>".dat","Table"];

Ppbc[order,temp]=ListPlot[data,Joined->True,PlotStyle->ColorData["Rainbow"][(Length[listtemp]-i-1)/(Length[listorder]-1)],PlotRange->All,ImageSize->700];

Ppbc2[order,temp]=ListPlot[data,Joined->True,PlotStyle->ColorData["Rainbow"][(Length[listtemp]-j-1)/(Length[listtemp]-1)],PlotRange->All,ImageSize->700],

{j,1,Length[listtemp]}],
{i,1,Length[listorder]}]

temp=0.0125;
order=18;
bla=Show[Ppbc[order,temp],Pobc[18,temp],ImageSize->800,PlotRange->All]