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]]];

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];