code for https://mathematica.stackexchange.com/q/193218

t = 1; nk = 1; J = 1; \[Lambda] = 0.3; \[Alpha] = 0.3;
Kx = -\[Pi]; Kx1 = \[Pi]; Ky = -\[Pi]; Ky1 = \[Pi];
sp1 = 0.01; d = 30;
G[A0_, \[Omega]0_] := (
  \[Omega] = \[Omega]0;
  T = 2 Pi/\[Omega];
  ClearAll[H, Hk, F];
  A = A0;
  E1[l_] := (\[Phi] = l \[Pi]; {Cos[\[Phi]], Sin[\[Phi]]}/2);
  E2[l_] := (\[Phi] = l \[Pi] + \[Pi]/2; {Cos[\[Phi]], Sin[\[Phi]]}/2);
  E\[Alpha][
    l_] := (\[Phi] =
     l \[Pi]/2 + \[Pi]/4; {Cos[\[Phi]], Sin[\[Phi]]} Sqrt[2]/2);
  H[k_] := H[k] = Sum[({
        {0,
         J BesselJ[k, A/2] Exp[-I {kx, ky}.E1[l]] Exp[
           I k (l \[Pi] + \[Pi]/2)],
         J BesselJ[k, A/2] Exp[-I {kx, ky}.E2[l]] Exp[I k (l \[Pi])]},
        {J BesselJ[k, -A/2] Exp[I {kx, ky}.E1[l]] Exp[
           I k (l \[Pi] + \[Pi]/2)], 0, 0},
        {J BesselJ[k, -A/2] Exp[I {kx, ky}.E2[l]] Exp[I k (l \[Pi])],
         0, 0}
       }), {l, 0, 1}] + Sum[({
        {0, \[Alpha] BesselJ[k, A/2] Exp[-I {kx, ky}.E1[l]] Exp[
           I k (l \[Pi] + \[Pi]/2)], \[Alpha] BesselJ[k,
           A/2] Exp[-I {kx, ky}.E2[l]] Exp[I k (l \[Pi])]},
        {\[Alpha] BesselJ[k, -A/2] Exp[I {kx, ky}.E1[l]] Exp[
           I k (l \[Pi] + \[Pi]/2)], 0, 0},
        {\[Alpha] BesselJ[k, -A/2] Exp[I {kx, ky}.E2[l]] Exp[
           I k (l \[Pi])], 0, 0}
       }), {l, 0}] - ({
       {0, \[Alpha] BesselJ[k, A/2] Exp[-I {kx, ky}.E1[1]] Exp[
          I k (1 \[Pi] + \[Pi]/2)], \[Alpha] BesselJ[k,
          A/2] Exp[-I {kx, ky}.E2[1]] Exp[I k (1 \[Pi])]},
       {\[Alpha] BesselJ[k, -A/2] Exp[I {kx, ky}.E1[1]] Exp[
          I k (1 \[Pi] + \[Pi]/2)], 0, 0},
       {\[Alpha] BesselJ[k, -A/2] Exp[I {kx, ky}.E2[1]] Exp[
          I k (1 \[Pi])], 0, 0}
      })+ Sum[({
        {0, 0, 0},
        {0,
         0, -I \[Lambda] BesselJ[k,
           Sqrt[2] A/2] Exp[-I {kx, ky}.E\[Alpha][l]] Exp[
           I k (l \[Pi]/2 - \[Pi]/4)]},
        {0,
         I \[Lambda] BesselJ[k, -Sqrt[2] A/2] Exp[
           I {kx, ky}.E\[Alpha][l]] Exp[I k (l \[Pi]/2 - \[Pi]/4)], 0}
       }), {l, 0, 3, 3}] - Sum[({
        {0, 0, 0},
        {0,
         0, -I \[Lambda] BesselJ[k,
           Sqrt[2] A/2] Exp[-I {kx, ky}.E\[Alpha][l]] Exp[
           I k (l \[Pi]/2 - \[Pi]/4)]},
        {0,
         I \[Lambda] BesselJ[k, -Sqrt[2] A/2] Exp[
           I {kx, ky}.E\[Alpha][l]] Exp[I k (l \[Pi]/2 - \[Pi]/4)], 0}
       }), {l, 1, 2}];

  Hk = Table[
    H[j - i] + IdentityMatrix[3] If[i == j, i \[Omega], 0], {i, -nk,
     nk}, {j, -nk, nk}];
  Hk = ArrayFlatten[Hk];
  F[x_, y_] := (
    SA = Eigensystem[N[Hk /. kx -> x /. ky -> y]];
    SA = Table[{Chop[SA[[1]][[i]]], SA[[2]][[i]]}, {i, (2 nk + 1) 3}];
    SA = Sort[SA, #1[[1]] < #2[[1]] &];
    SA1 = SA[[1 ;; (3 nk + 1)]][[All, 2]]);
  s1 = Table[
    F[x, y], {x, Kx + sp1, Kx1 + sp1, (Kx1 - Kx)/d}, {y, Ky + sp1,
     Ky1 + sp1, (Ky1 - Ky)/d}];
  s2 = Total[
     Table[Im[
       Log[Diagonal[(Conjugate[s1[[i, j]]].Transpose[
             s1[[i + 1, j]]]) (Conjugate[s1[[i + 1, j]]].Transpose[
             s1[[i + 1, j + 1]]]) (Conjugate[
             s1[[i + 1, j + 1]]].Transpose[
             s1[[i, j + 1]]]) (Conjugate[s1[[i, j + 1]]].Transpose[
             s1[[i, j]]])]]], {i, d}, {j, d}], 2]/(2 Pi);
  -Total[s2])
G[12, 10]