Clear[" \
Global`*"];

M = 4 ;   
v = {{1, 1, 0, 0}, {1, 1, 0, 0}};
FullForm[v]
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}])

v[[1]]

{1, 1, 0, 0}
index

FullForm[index]

??plus

FullForm[plus[_]]

FullForm[v]

v[[2]][[3]]

1
TreeForm[v]

1
F

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

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[factor2_ *s[arg1_], factor_. *s[arg2_]] := 
          factor2*factor*If[arg1 == arg2, 1, 0];

v
c[1, 2][s[v]]
scalarp[s[v], c[1, 2][s[v]]]

(* 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}]];
K[t_, 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[u_, vector_] := 
        Expand[          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*SS1[k][s[arg1]] + 
            factor2*SS1[k][s[arg2]]);
SS1[k_][
          factor1_. *s[arg1_]*factor2_. *
            s[arg2_]] := (factor1*SS1[k][s[arg1]]*
            factor2*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*SS2[k][s[arg1]] + 
            factor2*SS2[k][s[arg2]]);
SS2[k_][
          factor1_. *s[arg1_]*factor2_. *
            s[arg2_]] := (factor1*SS2[k][s[arg1]]*
            factor2*SS2[k][
              s[arg2]]);

SS2[k_][0] := 0;

SSz[k_][factor_. *s[arg_]] :=  
        factor*(Sz[k][Sz[k + 1][
                s[arg]]]);
SSz[k_][
          factor1_. *s[arg1_] + 
            factor2_. *s[arg2_]] := (factor1*SS2[k][s[arg1]] + 
            factor2*SS2[k][s[arg2]]);
SSz[k_][
          factor1_. *s[arg1_]*factor2_. *
            s[arg2_]] := (factor1*SS2[k][s[arg1]]*
            factor2*SS2[k][
              s[arg2]]);

SSz[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 } ]]]] ;
 gt[tt_, 
          uu_] := Sort[
            Thread[ Eigensystem[ 
                N[ h /.  {t -> tt,  
                      u  ->  uu } ]]]] [[1, 
            2]];
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}]); 
SMz[
          k_] := (slistz = Table[SSz[k][s[index[[j]]]], {j, MaxL}]; 
          Table[scalarp[s[index[[i]]], slistz[[j]]], {i, MaxL}, {j, 
              MaxL}]); 
KEM[t_] := (klist0 = Table[K[t, s[index[[j]]]], {j, MaxL}]; 
          Table[scalarp[s[index[[i]]], klist0[[j]]], {i, MaxL}, {j, 
              MaxL}]); 
UEM[u_] := (ulist0 = Table[U[u, s[index[[j]]]], {j, MaxL}]; 
          Table[scalarp[s[index[[i]]], ulist0[[j]]], {i, MaxL}, {j, 
              MaxL}]);

g[1. ]

SSp1[u_] := (GU = g[u]; 
          SSM1 = Sum[SM1[i], {i, M}];   Sum[ 
              GU[[j1]]  GU[[j2]] * SSM1[[j1, j2]], {j1, 
                MaxL}, {j2, MaxL}]/M ); 
SSp2[u_] := (GU = g[u]; 
          SSM2 = Sum[SM2[i], {i, M}];   Sum[ 
              GU[[j1]]  GU[[j2]] * SSM2[[j1, j2]], {j1, 
                MaxL}, {j2, MaxL}]/M );
\(\(SSpz[u_] := \((GU = g[u]; 
          SSMz = Sum[SMz[i], {i, M}]; \ \ Sum[\ 
              GU[\([j1]\)]\ \ GU[\([j2]\)]\ *\ SSMz[\([j1, j2]\)], {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;
KE[t_, u_] := (GU = gt[t, u]; kin =  KEM[t]; 
          Sum[ GU[[j1]]  GU[[j2]] * kin[[j1, j2]], {j1, 
                MaxL}, {j2, MaxL}]/M  ); 
UE[u_] := (GU = g[u];  pot = UEM[u];  
          Sum[ GU[[j1]]  GU[[j2]] * pot[[j1, j2]], {j1, 
                MaxL}, {j2, MaxL}]/M );
UEt[t_, u_] := (GU = gt[t, u];  pot = UEM[u];  
          Sum[ GU[[j1]]  GU[[j2]] * pot[[j1, j2]], {j1, 
                MaxL}, {j2, MaxL}]/M );

KE[1, 0] // N 
KE[1, 4] // N
KE[1, 8] // N
KE[1, 100] // N

UE[0] // N
UE[1] // N
UE[8] // N
UE[100] // N

SSp1[0] // N
SSp1[4] // N
SSp1[8] // N
SSp1[100] // N
SSpz[4]*3 // N
SSpz[8]*3 // N
SSpz[100]*3//N

SSp2[0] // N
SSp2[4] // N
SSp2[8] // N
SSp2[100] // N

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

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

Umax = 50;
Plot[dd[u], {u, 0, Umax}, 
      PlotPoints -> 6]
Plot[SSp1[u], {u, 0, Umax}, 
      PlotPoints -> 6]
Plot[SSp2[u], {u, 0, Umax}, PlotPoints -> 6]

1957  kg

Plot[KE[1, u], {u, 0, 50}, 
      PlotPoints -> 11]
Plot[UE[u], {u, 0, 50}, PlotPoints -> 11]