kp Hamiltonian generator (direct-decomposition method) 
Ref: arXiv:2104.08493

SetDirectory[NotebookDirectory[]];
Needs["MSGv8`"]
Names["MSGv8`*"]

Table[{g,PGName[[g]]},{g,PG//Length}]

MatNameD6h[[1;;64]]

Spinful 2-dim rep (inversion is represented by the identity matrix)

n=2;
{σ0,σx,σy,σz}=Table[PauliMatrix[i],{i,0,3}];
Drep=Table[MatrixExp[-I*kSO3RotAngle[[g]]*kSO3RotAxis[[g]].{σx,σy,σz}/2],{g,Uni}]//Simplify;
(*Basis of 2×2 Hermite matrices*)
X={σ0,σx,σy,σz};
G=9;
Gr=PG[[G]];
Print["Point group = ",PGName[[G]]];
Print["rep = ",Table[If[g<49,{MatNameOh[[g]],Drep[[g]]//MatrixForm},{MatNameD6h[[g]],Drep[[g]]//MatrixForm}],{g,Gr}]];
cs=Array[c,{n^2,n^2}];
Mrep=Table[Table[c[i,j]/.Solve[Drep[[g]].X[[j]].Inverse[Drep[[g]]]==Sum[X[[i]]*c[i,j],{i,n^2}],Flatten[cs[[All,j]]]]//First,{i,n^2},{j,n^2}],{g,Gr}]//Simplify;
Print["M = ",MatrixForm/@Mrep];
(*Basis of m-th order monomial *)
Do[
Print["---------- m = ",m," ----------"];
Nm=(m+2)*(m+1)/2;
mlis=Tuples[Table[i,{i,0,m}],3];
mlis=Select[mlis,Total[#]==m&];
f=Table[kx^(a[[1]])*ky^(a[[2]])*kz^(a[[3]]),{a,mlis}];
cs=Array[c,{Nm,Nm}];
Nrep=Table[Table[c[i,j]/.Solve[Table[D[(f[[j]]/.{kx->(kMat[[g]].{kx,ky,kz})[[1]],ky->(kMat[[g]].{kx,ky,kz})[[2]],kz->(kMat[[g]].{kx,ky,kz})[[3]]})-Sum[f[[i]]*c[i,j],{i,Nm}],{kx,a[[1]]},{ky,a[[2]]},{kz,a[[3]]}]==0,{a,mlis}],Flatten[cs[[All,j]]]]//First,{i,Nm},{j,Nm}],{g,Gr}]//Simplify;
Print["N = ",MatrixForm/@Nrep];
(*Intersection of nullspace*)
Lrep=Table[KroneckerProduct[Mrep[[g]],Nrep[[g]]],{g,Gr//Length}];
sols=NullSpace[Join@@Table[Lrep[[g]]-IdentityMatrix[n^2*Nm],{g,Gr//Length}]];
map[i_,a_]:=Nm*(i-1)+a;
Do[Print[Sum[sol[[map[i,a]]]*X[[i]]*f[[a]],{i,n^2},{a,Nm}]//MatrixForm],{sol,sols}];
,{m,{1,2,3}}];

Spinful 4-dim rep (inversion is represented by τz)

n=4;
{σ0,σx,σy,σz}=Table[KroneckerProduct[PauliMatrix[0],PauliMatrix[i]],{i,0,3}];
{τ0,τx,τy,τz}=Table[KroneckerProduct[PauliMatrix[i],PauliMatrix[0]],{i,0,3}];
Drep=Table[If[kO3Det[[g]]==1,τ0,τz].MatrixExp[-I*kSO3RotAngle[[g]]*kSO3RotAxis[[g]].{σx,σy,σz}/2],{g,Uni}]//Simplify;
(*Basis of 2×2 Hermite matrices*)
X=Flatten[Table[KroneckerProduct[PauliMatrix[i],PauliMatrix[j]],{i,0,3},{j,0,3}],1];
G=35;
Gr=PG[[G]];
Print["Point group = ",PGName[[G]]];
Print["rep = ",Table[If[g<49,{MatNameOh[[g]],Drep[[g]]//MatrixForm},{MatNameD6h[[g]],Drep[[g]]//MatrixForm}],{g,Gr}]];
cs=Array[c,{n^2,n^2}];
Mrep=Table[Table[c[i,j]/.Solve[Drep[[g]].X[[j]].Inverse[Drep[[g]]]==Sum[X[[i]]*c[i,j],{i,n^2}],Flatten[cs[[All,j]]]]//First,{i,n^2},{j,n^2}],{g,Gr}]//Simplify;
Print["M = ",MatrixForm/@Mrep];
(*Basis of m-th order monomial *)
Do[
Print["---------- m = ",m," ----------"];
Nm=(m+2)*(m+1)/2;
mlis=Tuples[Table[i,{i,0,m}],3];
mlis=Select[mlis,Total[#]==m&];
f=Table[kx^(a[[1]])*ky^(a[[2]])*kz^(a[[3]]),{a,mlis}];
cs=Array[c,{Nm,Nm}];
Nrep=Table[Table[c[i,j]/.Solve[Table[D[(f[[j]]/.{kx->(kMat[[g]].{kx,ky,kz})[[1]],ky->(kMat[[g]].{kx,ky,kz})[[2]],kz->(kMat[[g]].{kx,ky,kz})[[3]]})-Sum[f[[i]]*c[i,j],{i,Nm}],{kx,a[[1]]},{ky,a[[2]]},{kz,a[[3]]}]==0,{a,mlis}],Flatten[cs[[All,j]]]]//First,{i,Nm},{j,Nm}],{g,Gr}]//Simplify;
Print["N = ",MatrixForm/@Nrep];
(*Intersection of nullspace*)
Lrep=Table[KroneckerProduct[Mrep[[g]],Nrep[[g]]],{g,Gr//Length}];
sols=NullSpace[Join@@Table[Lrep[[g]]-IdentityMatrix[n^2*Nm],{g,Gr//Length}]];
map[i_,a_]:=Nm*(i-1)+a;
Do[Print[Sum[sol[[map[i,a]]]*X[[i]]*f[[a]],{i,n^2},{a,Nm}]//MatrixForm],{sol,sols}];
,{m,{1,2,3}}];