Clear[RawHamBloch,liens,cparmodcom];

Octonacci[k_]:=Module[{octo,i},
octo=Table[0,{i,1,k}];
octo[[1]]=1;
octo[[2]]=2;
Do[ octo[[i]]=octo[[i-2]]+2*octo[[i-1]],{i,3,k}];
octo];
locto=Octonacci[12];

locto

liens=Import["liensPBC8.dat"];
cparmodcom=Import["PtsPBC8.dat"];
l=Length[liens];
ll=Length[cparmodcom];

cut=0.5;


RawHamBloch=Reap[Do[{rat=N[locto[[9]]*Sqrt[2]];dx=cparmodcom[[liens[[m,1]],1]]-cparmodcom[[liens[[m,2]],1]];dy=cparmodcom[[liens[[m,1]],2]]-cparmodcom[[liens[[m,2]],2]];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dx]>cut)&&(Abs[dy]>cut),(*Print["Cond 1 "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
Sow[{{liens[[m,1]],liens[[m,2]]},- Exp[I*rat*(kx*(dx-Sign[dx])+ky*(dy-Sign[dy]))]}];Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-I*rat*(kx*(dx-Sign[dx])+ky*(dy-Sign[dy]))]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dx]>cut)&&(Abs[dy]<cut),(*Print["Cond 2 "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[I*rat*(kx*(dx-Sign[dx])+ky*(dy))]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-I*rat*(kx*(dx-Sign[dx])+ky*(dy))]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dy]>cut)&&(Abs[dx]<cut),(*Print["Cond 3 "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[I*rat*(kx*(dx)+ky*(dy-Sign[dy]))]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-I*rat*(kx*(dx)+ky*(dy-Sign[dy]))]}];];

If[(liens[[m,1]]<liens[[m,2]])&&(Abs[dy]<cut)&&(Abs[dx]<cut),(*Print["Cond 4 "," ",dx," ",dy," ",liens[[m,1]]," ",liens[[m,2]]];*)
Sow[{{liens[[m,1]],liens[[m,2]]},-Exp[I*rat*(kx*(dx)+ky*(dy))]}];
Sow[{{liens[[m,2]],liens[[m,1]]},-Exp[-I*rat*(kx*(dx)+ky*(dy))]}];]
},{m,1,l}]][[2,1]];

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

sp9xv2=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*N[Pi/(8*985*Sqrt[2])],ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,20,Method->{"Arnoldi",MaxIterations->4000}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];

Export["sp9x.dat",N[sp9xv2]];

Length[sp9xv2]

sp9yv2=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={ky->(n-3)*N[Pi/(8*985*Sqrt[2])],kx->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,20,Method->{"Arnoldi",MaxIterations->4000}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];
sp9xyv2=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*N[Pi/(8*985*Sqrt[2])],ky->(n-3)*N[Pi/(8*985*Sqrt[2])]};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,20,Method->{"Arnoldi",MaxIterations->4000}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];
Export["sp9y.dat",N[sp9yv2]];
Export["sp9xy.dat",N[sp9xyv2]];

Show[ListPlot[Transpose[{Transpose[sp9xv2][[1]],Transpose[sp9xv2][[2]]}],PlotRange->{{-0.0006,0.0006},{-4.22170,-4.22169}}],Plot[-4.221696440289577`+1.4*(x^2),{x,-0.0006,0.0006}]]

1

Fit[Take[Transpose[{Transpose[sp9xv2][[1]],Transpose[sp9xv2][[2]]}],{2,4}],{1,x,x^2,x^3},x]

Show[ListPlot[Transpose[{Transpose[sp9yv2][[1]],Transpose[sp9yv2][[2]]}],PlotRange->{{-0.0006,0.0006},{-4.22170,-4.22169}}],Plot[-4.221696440289577`+1.4*(x^2),{x,-0.0006,0.0006}]]

Fit[Take[Transpose[{Transpose[sp9yv2][[1]],Transpose[sp9yv2][[2]]}],{2,4}],{1,x,x^2,x^3},x]

Fit[Take[Transpose[{Transpose[sp9xyv2][[1]],Transpose[sp9xyv2][[2]]}],{2,4}],{1,x,x^2,x^3},x]

rulef={ky->2*N[Pi/(8*985*Sqrt[2])],kx->0};
ham=SparseArray[RawHamBloch/.rulef];

eig8000=Eigenvalues[ham,20,Method->{"Arnoldi",MaxIterations->8000}];

i fe

eig8000

{-4.221696959957663` ,4.221686390481906`,4.221679793170643`,4.221676138708123`,-4.221670331066098`,-4.22166194560218`,-4.2216448656372005`,4.221636838833769`,-4.221619435632939`,4.221615100282034`,-4.221571829336112`,4.221562776382761`,-4.221534683598248`,4.221521041935966`,4.221466199678775`,-4.221452542731292`,-4.221384379564471`,4.221378911397559`,4.221290546711841`,-4.221259154976075`}

Transpose[{Transpose[sp9yv2][[1]],Transpose[sp9yv2][[2]]}]

{{-0.0005638176317459855`,-4.221694945376788`},{-0.00028190881587299274`,-4.221695968307904`},{0.`,-4.2216964931114225`},{0.00028190881587299274`,-4.221696485160215`},{0.0005638176317459855`,-4.221695798209291` }}

sp9xv3=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*N[Pi/(8*985*Sqrt[2])],ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,20,Method->{"Arnoldi",MaxIterations->800}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];

sp9yv3=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={ky->(n-3)*N[Pi/(8*985*Sqrt[2])],kx->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,20,Method->{"Arnoldi",MaxIterations->8000}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];

sp9xyv3=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*N[Pi/(8*985*Sqrt[2])],ky->(n-3)*N[Pi/(8*985*Sqrt[2])]};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,20,Method->{"Arnoldi",MaxIterations->4000}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];

Export["sp9yv3.dat",N[sp9yv3]];

sp9xv4=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*N[Pi/(8*985*Sqrt[2])],ky->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,40,Method->{"Arnoldi",MaxIterations->10000}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];


sp9yv4=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={ky->(n-3)*N[Pi/(8*985*Sqrt[2])],kx->0};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,40,Method->{"Arnoldi",MaxIterations->10000}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];

sp9xyv4=Monitor[Reap[
For[n=1,n≤5,n++,
rulef={kx->(n-3)*N[Pi/(8*985*Sqrt[2])],ky->(n-3)*N[Pi/(8*985*Sqrt[2])]};
ham=SparseArray[RawHamBloch/.rulef];
Sow[Flatten[{(n-3)*N[Pi/(985*8*Sqrt[2])],Sort[Eigenvalues[ham,40,Method->{"Arnoldi",MaxIterations->10000}]]}]/.rulef]
(*,Method->{"Arnoldi","MaxIterations"->4000}]
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]],n];

Dimensions[sp9xv4[[2]]]

Show[ListPlot[Transpose[{Transpose[sp9yv4][[1]],Transpose[sp9yv4][[2]]}],PlotRange->{{-0.0006,0.0006},{-4.221701,-4.22169}}]]

sp9yv3=Import["sp9yv3.dat"];

Fit[Transpose[{Transpose[sp9xv3][[1]],Transpose[sp9xv3][[2]]}],{1,x,x^2,x^3},x]

-4.22169741477318`+3.3923846007425438`*^-6 x+1.4293031615136789`  (x^2)+2.9750366804075763` (x^3)

Fit[Transpose[{Transpose[sp9xv3][[1]],Transpose[sp9xv3][[2]]}],{1,x^2},x]

-4.221697414773179`+1.4293031665871567` (x^2)

-4.221697414773179`+1.4293031665871567` ((x )^2)+14.139126788118853` (x^3)

Fit[Transpose[{Transpose[sp9yv3][[1]],Transpose[sp9yv3][[2]]}],{1,x,x^2,x^3},x]

Fit[Transpose[{Transpose[sp9xyv3][[1]],Transpose[sp9xyv3][[2]]}],{1,x,x^2,x^3},x]

Fit[Transpose[{Transpose[sp9xyv3][[1]],Transpose[sp9xyv3][[2]]}],{1,x^2},x]

alpha8x=2*1.4293031615136789`
alpha8y=2*1.4433801914743851`
alpha8xy=2.829256858915674`-1.4293031615136789`-1.4433801914743851`

2.8586063230273577`

2.8867603829487702`

-0.043426494072389854`

minv8=Sqrt[alpha8x*alpha8y-alpha8xy*alpha8xy]

2.872320598988701`

minv8v2=Sqrt[4*1.4293031665871567`*1.4433801973496132`-(2.8292568648934093`-1.4293031665871567`-1.4433801973496132`)^2]

Fit[Transpose[{Transpose[sp9xv4][[1]],Transpose[sp9xv4][[2]]}],{1,x,x^2,x^3},x]
Fit[Transpose[{Transpose[sp9yv4][[1]],Transpose[sp9yv4][[2]]}],{1,x,x^2,x^3},x]
Fit[Transpose[{Transpose[sp9xyv4][[1]],Transpose[sp9xyv4][[2]]}],{1,x,x^2,x^3},x]

alpha8x4=2*1.394566257997868`
alpha8y4=2*1.3945622204518224`
alpha8xy4=2.7879823947115145`-1.394566257997868`-1.3945622204518224`

minv8v4=Sqrt[alpha8x4*alpha8y4-alpha8xy4*alpha8xy4]