Clear[" Global`*"];

M = 4 ;   
v = {{1, 1, 0, 0}, {1, 1, 0, 0}};
MP = 
        Length[Table[Permutations[v[[1]]]]];
index = 
        Flatten[Table[{Permutations[v[[1]]][[i]], Permutations[
                  v[[2]]][[j]]}, {i, MP}, {j, MP}], 
          1];
plus[k_, sigma_][arg_] := 
        ReplacePart[arg, 1, {sigma, k}];
minus[k_, sigma_][arg_] := 
        ReplacePart[arg, 0, {sigma, k}];
signn[1, 1][arg_] := 1
signn[k_, 1][arg_] := (-1)^
        Sum[arg[[1, i]], {i, 1, k - 1}]
signn[1, 2][arg_] :=   (-1)^
        Sum[arg[[1, i]], {i, 1, M}]
signn[k_, 2][
        arg_] := (-1)^(Sum[arg[[2, i]], {i, 1, k - 1}] + 
            Sum[arg[[1, i]], {i, 1, M}])

c[k_, sigma_][0] := 0;
c[k_, sigma_][
        factor_. *s[arg_]] := 
      factor*(arg[[sigma, k]])*signn[k, sigma][arg]*
        s[minus[k, sigma][arg]]
cdg[k_, sigma_][0] := 0;
cdg[k_, sigma_][factor_. *s[arg_]] := 
          factor*(1 - arg[[sigma, k]])*signn[k, sigma][arg]*
            s[plus[k, sigma][arg]];
(* cdg1[k_, sigma_][0] := 0;
cdg1[k_, sigma_][factor_. *arg_] := 
            factor*(1 - arg[[sigma, k]])*signn[k, sigma][
                arg]*plus[k, sigma][arg] *)
c[M + 1, sigma_] := 
      c[1, sigma]
cdg[M + 1, sigma_] := cdg[1, sigma]

c[M + 2, sigma_] := 
      c[2, sigma]
cdg[M + 2, sigma_] := cdg[2, sigma]

n[k_, sigma_][0] := 
      0
n[k_, sigma_][factor_. *s[arg_]] := 
        factor*arg[[sigma, k]]*s[arg]

scalarp[a_, 0] := 
        0;
scalarp[a_, b_ + c_] := 
        scalarp[a, b] + scalarp[a, c];
scalarp[
          s[arg1_], factor_. *s[arg2_]] := 
        factor*If[arg1 == arg2, 1, 0];

(* 
      The Hamiltonian *) 
H[vector_] := 
        Expand[            -t*
              Sum[cdg[k, sigma][ c[k + 1, sigma][
                      vector] ]  + cdg[k + 1, sigma][ c[k, sigma][
                      vector] ], 
                                       {k, M}, {sigma, 
                  2}] + 
                                       u*
              Sum[n[k, 1][n[k, 2][vector]], {k, 1, 
                  M}]];

Sp[k_][
          factor_. *
            s[arg_]] :=           factor*(cdg[k, 
                1][ c[k, 2][s[arg]]]);
Sm[k_][
          factor_. *
            s[arg_]] :=          factor*(cdg[k, 
                2][ c[k, 1][s[arg]]]);
Sz[k_][
          factor_. *
            s[arg_]] :=          factor*(cdg[k, 
                    1][ c[k, 1][s[arg]]] - cdg[k, 2][ c[k, 2][
                    s[arg]]])/2;
Sp[k_][0] := 0;
Sm[k_][0] := 0;
Sz[k_][0] := 0;

SS1[k_][factor_. *s[arg_]] :=  
        factor*(0.5*(Sp[k][Sm[k + 1][s[arg]]] + Sm[
                        k][Sp[k + 1][s[arg]]]) + Sz[
                  k][Sz[k + 1][
                  s[arg]]]);
SS1[k_][factor1_. *s[arg1_] + factor2_. *s[arg2_]] := 
        factor1*factor2*(SS1[k][s[arg1]] + SS1[k][
                s[arg2]]);
SS1[k_][0] := 0;

SS2[k_][factor_. *s[arg_]] :=  
        factor*(0.5*(Sp[k][Sm[k + 2][s[arg]]] + Sm[
                        k][Sp[k + 2][s[arg]]]) + Sz[
                  k][Sz[k + 2][
                  s[arg]]]);
SS2[k_][factor1_. *s[arg1_] + factor2_. *s[arg2_]] := 
        factor1*factor2*(SS2[k][s[arg1]] + SS2[k][
                s[arg2]]);
SS2[k_][0] := 0;

MaxL = Length[index];
h = (hlist = Table[H[s[index[ [j]] ]], {j, MaxL}]; 
          Table[scalarp[  s[ index[[i]] ], hlist[[j]] ],  {i, 
              MaxL}, {j, MaxL}]);
  g[
          uu_] := Sort[
            Thread[ Eigensystem[ 
                N[ h /.  {t -> 1.0,  
                      u  ->  uu } ]]]] [[1, 
            2]];
gs[uu_] := 
        Sort[Thread[ 
            Eigensystem[ 
              N[ h /.  {t -> 1.0,  
                    u  ->  uu } ]]]] ;
SM1[
          k_] := (slist1 = Table[SS1[k][s[index[[j]]]], {j, MaxL}]; 
          Table[scalarp[s[index[[i]]], slist1[[j]]], {i, MaxL}, {j, 
              MaxL}]);
SM2[
          k_] := (slist2 = Table[SS2[k][s[index[[j]]]], {j, MaxL}]; 
          Table[scalarp[s[index[[i]]], slist2[[j]]], {i, MaxL}, {j, 
              MaxL}]);

\(SSp1[u_] := 
      Sum[\ \(g[u]\)[\([j1]\)]\ \ \(g[u]\)[\([j2]\)]\ *\ 
            \(SM1[i]\)[\([j1, j2]\)], {i, M}, {j1, MaxL}, {j2, MaxL}]/M; \ 

    


    SSp2[u_] := 
      Sum[\ \(g[u]\)[\([j1]\)]\ \ \(g[u]\)[\([j2]\)]\ *\ 
            \(SM2[i]\)[\([j1, j2]\)], {i, M}, {j1, MaxL}, {j2, MaxL}]/M; \ 

    


    dd[u_] := 
      Sum[Sum[Abs[\(g[u]\)[\([j]\)]]^2*index[\([j, 1, k]\)]*
              index[\([j, 2, k]\)], {j, MaxL}], {k, M}]/M; \)

MatrixForm[h]
MatrixForm[gs[0]]
g[0]
Eigenvalues[h]

Umax = 
      6. 
Plot[dd[u], {u, 0, Umax}]
Plot[SSp1[u], {u, 0, Umax}]
Plot[SSp2[u], {u, 0, Umax}]

(* =============================================================== *) 
TeXForm[Eigenvalues[h]] >> 
        hubbard.tex;