Continuum model Hamiltonian for TMDS

Question

I am trying to implement the continuum model Hamiltonian for transition metal dichalcogenide systems as described in the reference: PRL 122, 086402 (2019)

   (* Define the constants*)
   (* twist angle*)
   \[Theta] = 1.2 Degree;
   (* For MoTe2 V, psi, w*)
   V = 8 (* meV*);
   psi = -89.6 ;
   w = -8.5 (*meV*);
   Vz = -2 (* displacement field meV*);
   meff = 0.62 (* valence band effective mass*);
   a0 = 3.472 (*angstrom*);
   NN = 6(*truncate range*);

   \[Theta] = 
   1.2 Degree; (*\[Theta]\[Pi]/180conversion from degress to radian*)

    V = V/1000;
    w = w/1000;
    Vz = Vz/1000;
    psi = psi Degree;
    hbar = 4.135667*10^-15/2/\[Pi];
    me = 0.51099895*10^6;
    meff = meff*me;
    kin = -hbar^2/2/meff*9*10^16;
    a0 = 3.472 ;
    a0 = a0*10^-10;
    a1 = a0 {1, 0};
    a2 = a0 {Cos[\[Pi]/3], Sin[\[Pi]/3]};
    (* Reciprocal lattice vectors*)
    G1 = (4 \[Pi])/((Sqrt[3]) a0) {0, 1};
    G2 = (4 \[Pi])/((Sqrt[3]) a0) {-Sin[\[Pi]/3], Cos[\[Pi]/3]};
    (*G2=(4\[Pi])/((Sqrt[3])a0){Cos[\[Pi]/3],Sin[\[Pi]/3]};*)

    g1 = \[Theta] {G1[[2]], -G1[[1]]};
    g2 = \[Theta] {G2[[2]], -G2[[1]]};
    kD = g1[[1]]/Sqrt[3];

    Kb = {-kD*Sqrt[3]/2, -kD/2}; (* K botom layer *)
    Kt = {-kD*Sqrt[3]/2, 
     kD/2} ;   (* K top layer *)

     (* Define the lattice*)
      n = 6;
      nlayer = (2 n + 1)^2;
      L = Table[0, {i, 1, 2 (2 n + 1)^2}, {j, 1, 2}];
     invL = Table[0, {i, 1, 2 n + 1}, {j, 1, 2 n + 1}];
     k = 1;
     count = 0;
     Do[Do[
    (*bottom layer*)
     L[[k, 1]] = i; L[[k, 2]] = j;
     invL[[i + n + 1, j + n + 1]] = count;
     count = count + 1;

     (* top layer*)
     L[[k + nlayer, 1]] = i; L[[k + nlayer, 2]] = j;
     k = k + 1;
     , {j, -n, n}], {i, -n, n}]

      (* check the form of inVL*)
       MatrixForm[invL];

       (* Hamiltonian implementation *)
        Clear[kx, ky]

         HHHH[kx_, ky_] :=
          Module[{ h, NN},

      NN = 6;
      siteN = (2 NN + 1)^2;
      h = ConstantArray[0, {2 siteN, 2 siteN}];
     Do[ 
      ix = L[[i, 1]];
      iy = L[[i, 2]];

      (* diagnal elements in the Hamiltonian*)
  
       ax = kx - Kb[[1]] + ix g1[[1]] + iy g2[[1]];
         ay = ky - Kb[[2]] + ix g1[[2]] + iy g2[[2]];
         h[[i, i]] += kin (ax^2 + ay^2) + Vz/2;

       ax = kx - Kt[[1]] + ix g1[[1]] + iy g2[[1]];
         ay = ky - Kt[[2]] + ix g1[[2]] + iy g2[[2]];
         h[[i + siteN, i + siteN]] += kin (ax^2 + ay^2) - Vz/2;

        h[[i + siteN, i]] += w;
         h[[i, i + siteN]] += w;

        If[ix != NN,
         j = invL[[ix + 1 + NN + 1, iy + NN + 1]];
          h[[j + 1, i]] += V Exp[I*psi];
          h[[j + siteN + 1, i + siteN]] += V Exp[-I*psi];
         ];

          If [iy != -NN,
          j = invL[[ix + NN + 1, iy - 1 + NN + 1]];
                 h[[j + 1, i]] += V Exp[I psi];
                 h[[j + siteN + 1, i + siteN]] += V*Exp[-I psi];
                 h[[j + 1, i + siteN]] += w;
           ];

             If [ix != -NN && iy != NN,

               j = invL[[ix - 1 + NN + 1, iy + 1 + NN + 1]];
                 h[[j + 1, i]] += V*Exp[I*psi];
                 h[[j + siteN + 1, i + siteN]] += V*Exp[-I*psi];
                 h[[j + siteN + 1, i]] += w;
                 ];

            If [ix != -NN,
             j = invL[[ix - 1 + NN + 1, iy + NN + 1]];
                 h[[j + 1, i]] += V*Exp[-I*psi];
                h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
              ];
          If [iy != NN,
          j = invL[[ix + NN + 1, iy + 1 + NN + 1]];
                 h[[j + 1, i]] += V*Exp[-I*psi];
                h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
                h[[j + siteN + 1, i]] += w;
              ];

         If [ix != NN && iy != -NN,

                 j = invL[[ix + 1 + NN + 1, iy - 1 + NN + 1]];
                 h[[j + 1, i]] += V*Exp[-I*psi];
                 h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
                 h[[j + 1, i + siteN]] += w;];

               , {i, 1, siteN}];

             Return[h]

              ]

                Num = 50;
                KXX1 = N@Subdivide[g1[[1]]/2, 0, Num - 1];
                 KXX2 = N@Subdivide[0, -g1[[1]]/2, Num - 1];
                KXX3 = N@Subdivide[-g1[[1]]/2, -g1[[1]]/2, Num - 1];
                KXX4 = N@Subdivide[-g1[[1]]/2, g1[[1]]/2, (2 Num - 1)];
                KX = Join[KXX1, KXX2, KXX3, KXX4];

                KYY1 = N@Subdivide[-g1[[1]]/(Sqrt [3]/2), 0, Num - 1];
               KYY2 = N@Subdivide[0, g1[[1]]/(Sqrt[3]/2), Num - 1];
                KYY3 = N@
                Subdivide[g1[[1]]/(Sqrt[3]/2), -g1[[1]]/(Sqrt [3]/2), Num - 1];
              KYY4 = N@
              Subdivide[-g1[[1]]/(Sqrt [3]/2), -g1[[1]]/(Sqrt [3]/2), (2 Num - 
                    1)];
              KY = Join[KYY1, KYY2, KYY3, KYY4];

           Eigen = Table[0, {i1, Length[KX]}, {j1, 2 siteN}];

          Do[Eigen[[k]] = Eigenvalues[HHHH[KX[[k]], KY[[k]]]] // Sort, {k, 1, 
          Length[KX]}];

        out = ParallelTable[
         Eigenvalues[HHHH[KX[[k]], KY[[k]]]] // Sort, {k, 1, Length[KX]}];
           ListLinePlot[Transpose[out], PlotRange -> {{0, 250}, {0, 0.025}}]

The resulting band structure is look like this

I recently came across the github code which reproduce the band structure from the reference. (https://github.com/zihaophys/twisted_bilayer_graphene/blob/master/TMD/homobilayer/moireTMD.py)

I have studied and trying to reproduce in exact same way. Can someone please point out how to fix the band structure.

Thanks