过程版

ClearAll["`*"]
sd1k={kx->kx/2-√3 *ky/2,ky->-(Sqrt[3])kx/2-ky/2,kz-> kz};C3zk={kx->-kx/2+√3 *ky/2,ky->-(Sqrt[3])kx/2-ky/2, kz->kz};shk={kx->kx,ky->ky,kz->-kz};
PTk={kx->kx,ky->ky,kz->kz};
sd1=-I*KroneckerProduct[PauliMatrix[0], PauliMatrix[1]];C3z=(-KroneckerProduct[PauliMatrix[0], PauliMatrix[0]]-I*Sqrt[3]*KroneckerProduct[PauliMatrix[3], PauliMatrix[3]])/2;sh=KroneckerProduct[PauliMatrix[0], PauliMatrix[3]];
PT=KroneckerProduct[PauliMatrix[1], PauliMatrix[0]];
UOperationRes={{C3z,C3zk},{sh,shk},{sd1,sd1k}};(*整合所有的幺正操作*)
AUOperationRes={{PT,PTk}};(*整合所有的反幺正操作*)
OperationRes={UOperationRes,AUOperationRes};(*将所有对称操作整合起来,总体上分为两类,幺正操作和非幺正操作*)
(*----------------------------------------------------------------------------------------*)
nk=1;(*哈密顿量多项式部分展开的阶数*)
Vk={kx,ky,kz};
varlist=MonomialList[Expand[(((1+Total[Vk]))^nk)]/._Integer*b_->b,"NegativeDegreeReverseLexicographic"];(*展开多项式并只提取多项式部分，扔掉系数*)
matdim=4;(*想要展开矩阵的维度*)
numvar=Length[varlist];(*哈密顿量展开,多项式的数量*)
tdc=Array[dc,{(Sum[i,{i,matdim}]+Sum[i,{i,matdim-1}])*numvar}];
matcon1[numvar_,matdim_]:=Module[{mdr,mdi},mdr=Block[{ta,im},ta=Table[IdentityMatrix[matdim],{j,1,numvar}];(*根据表示的维度构建单位矩阵，数量与多项式的数目相同*)im=0;
Do[Do[im=im+1;
ta[[i,ii,ij]]=dc[im],{ii,1,matdim},{ij,ii,matdim}],{i,1,numvar}];(*i索引ta保存的每一个单位矩阵，ii和ij则是单位矩阵中元素的索引*)Do[Do[ta[[j,ij,ii]]=ta[[j,ii,ij]],{ii,1,matdim},{ij,ii,matdim}],{j,1,numvar}];(*这里是对矩阵中的元素做出一些限制,也就是令ij位置元素和ji位置元素相同，即对称矩阵要求ij=ji*)ta];
mdi=Block[{ta,im},ta=Table[IdentityMatrix[matdim],{j,1,numvar}];(*这里根据多项式展开的数量来确定需要构建的单位矩阵的数量*)im=Sum[i,{i,matdim}]*numvar;(*这里就是构建一个对称矩阵和反对称矩阵来表示厄米的哈密顿量，而自然这两种矩阵中的参数是要不相同的，所以这里就是用im这个索引来区分这些参数，因此在这里的时候使用到的所以就是在上面一个矩阵索引的基础上继续进行叠加，所以这里也就是为什么tdc要构建的dc的数量是(Sum[i,{i,id}]+Sum[i,{i,id-1}])*imd个*)Do[Do[If[ii==ij,ta[[i,ii,ij]]=0,im=im+1;
ta[[i,ii,ij]]=I*dc[im]],{ii,1,matdim},{ij,ii,matdim}],{i,1,numvar}];(*此时限制对角元素为零，非对角和mdr相同只不过乘以了虚数I*)Do[Do[ta[[j,ij,ii]]=-ta[[j,ii,ij]],{ii,1,matdim},{ij,ii,matdim}],{j,1,numvar}];(*此时给出的是限制矩阵为反对称矩阵:ij=-ji*)ta];

{mdr,mdi}](*这里在构建对称矩阵和反对称矩阵的时候使用了一个参量dc来表示矩阵中的参数,用的是全局变量,因为后面要对这个参数进行限制*)
matcon[numvar_,matdim_]:=Module[{mdr,mdi},mdr=Block[{ta,im},ta=Table[IdentityMatrix[matdim],{j,1,numvar}];(*根据表示的维度构建单位矩阵，数量与多项式的数目相同*)im=0;
Do[Do[im=im+1;
ta[[i,ii,ij]]=tdc[[im]],{ii,1,matdim},{ij,ii,matdim}],{i,1,numvar}];(*i索引ta保存的每一个单位矩阵，ii和ij则是单位矩阵中元素的索引*)Do[Do[ta[[j,ij,ii]]=ta[[j,ii,ij]],{ii,1,matdim},{ij,ii,matdim}],{j,1,numvar}];(*这里是对矩阵中的元素做出一些限制,也就是令ij位置元素和ji位置元素相同，即对称矩阵要求ij=ji*)ta];
mdi=Block[{ta,im},ta=Table[IdentityMatrix[matdim],{j,1,numvar}];(*这里根据多项式展开的数量来确定需要构建的单位矩阵的数量*)im=Sum[i,{i,matdim}]*numvar;(*这里就是构建一个对称矩阵和反对称矩阵来表示厄米的哈密顿量，而自然这两种矩阵中的参数是要不相同的，所以这里就是用im这个索引来区分这些参数，因此在这里的时候使用到的所以就是在上面一个矩阵索引的基础上继续进行叠加，所以这里也就是为什么tdc要构建的dc的数量是(Sum[i,{i,id}]+Sum[i,{i,id-1}])*imd个*)Do[Do[If[ii==ij,ta[[i,ii,ij]]=0,im=im+1;
ta[[i,ii,ij]]=I*tdc[[im]]],{ii,1,matdim},{ij,ii,matdim}],{i,1,numvar}];(*此时限制对角元素为零，非对角和mdr相同只不过乘以了虚数I*)Do[Do[ta[[j,ij,ii]]=-ta[[j,ii,ij]],{ii,1,matdim},{ij,ii,matdim}],{j,1,numvar}];(*此时给出的是限制矩阵为反对称矩阵:ij=-ji*)ta];{mdr,mdi}]
hammat[numvar_,matdim_]:=Module[{mat1},mat1=matcon[numvar,matdim];
varlist.(mat1[[1]]+mat1[[2]])](*对称矩阵*)
hammatC[numvar_,matdim_]:=Module[{mat1},mat1=matcon[numvar,matdim];
varlist.(mat1[[1]]-mat1[[2]])](*反对称矩阵构建*)
diffham[opmat_,opk_,numvar_,matdim_]:=opmat.hammat[numvar,matdim].Inverse[opmat]-(hammat[numvar,matdim]/.opk)
diffhamC[opmat_,opk_,numvar_,matdim_]:=opmat.hammatC[numvar,matdim].Inverse[opmat]-(hammat[numvar,matdim]/.opk)
diffcoff[opmat_,opk_,numvar_,matdim_]:=Flatten[CoefficientList[diffham[opmat,opk,numvar,matdim],Vk]]
diffcoffC[opmat_,opk_,numvar_,matdim_]:=Flatten[CoefficientList[diffhamC[opmat,opk,numvar,matdim],Vk]]
constrain[opmat_,opk_,numvar_,matdim_]:=Module[{caplist(*这个参数我是想作为方程求解的时候的变量,因为通过对称性限制之后,一些参数可能会被限制为零,那么在方程求解过程中这些参量就不会作为变量来求解了*),opcons(*这个参量存储对称限制后相差的部分*),re1(*这个变量就用来存储系数方程的解*)},
opcons=diffcoff[opmat,opk,numvar,matdim];
re1=Flatten[Quiet@Solve[opcons==Table[0,{i,1,Length[opcons]}]]];tdc=(tdc/.re1)]


constrainC[opmat_,opk_,numvar_,matdim_]:=Module[{caplist(*这个参数我是想作为方程求解的时候的变量,因为通过对称性限制之后,一些参数可能会被限制为零,那么在方程求解过程中这些参量就不会作为变量来求解了*),opcons(*这个参量存储对称限制后相差的部分*),re1(*这个变量就用来存储系数方程的解*)},
opcons=diffcoffC[opmat,opk,numvar,matdim];
re1=Flatten[Quiet@Solve[opcons==Table[0,{i,1,Length[opcons]}]]];tdc=(tdc/.re1)]

constrainU[OpRes_,numvar_,matdim_]:=Module[{coff},Do[coff=constrain[OpRes[[i,1]],OpRes[[i,2]],numvar,matdim],{i,1,Length[OpRes]}];
coff](*遍历所有的幺正操作操作进行约束*)


constrainAU[OpRes_,numvar_,matdim_]:=Module[{coff},Do[coff=constrainC[OpRes[[i,1]],OpRes[[i,2]],numvar,matdim],{i,1,Length[OpRes]}];
coff](*遍历所有的反幺正操作操作进行约束*)


constrainAll[OpRes_,numvar_,matdim_]:=Module[{},constrainU[OpRes[[1]],numvar,matdim];constrainAU[OpRes[[2]],numvar,matdim];](*整合幺正与非幺正操作*)


constrainAll[OperationRes,numvar,matdim]
tdc=Flatten[tdc/.Quiet[Solve[Variables[tdc]==Table[(Subscript[c, i]),{i,1,Length[Variables[tdc]]}],Variables[tdc]]]];(*这里看起来是个方程求解,其实就是想将之前的构建矩阵时候的参量换个形式表示出来*)
Style[MatrixForm[FullSimplify[hammat[numvar,matdim]/.{kx->(Subscript[k, x]),ky->(Subscript[k, y]),kz->(Subscript[k, z])}]],FontFamily->"Times",20,Blue]

函数版

ClearAll["`*"]
sd1k={kx->kx/2-√3 *ky/2,ky->-(Sqrt[3])kx/2-ky/2,kz-> kz};C3zk={kx->-kx/2+√3 *ky/2,ky->-(Sqrt[3])kx/2-ky/2, kz->kz};shk={kx->kx,ky->ky,kz->-kz};
PTk={kx->kx,ky->ky,kz->kz};
sd1=-I*KroneckerProduct[PauliMatrix[0], PauliMatrix[1]];C3z=(-KroneckerProduct[PauliMatrix[0], PauliMatrix[0]]-I*Sqrt[3]*KroneckerProduct[PauliMatrix[3], PauliMatrix[3]])/2;sh=KroneckerProduct[PauliMatrix[0], PauliMatrix[3]];
PT=KroneckerProduct[PauliMatrix[1], PauliMatrix[0]];
UOperationRes={{C3z,C3zk},{sh,shk},{sd1,sd1k}};(*整合所有的幺正操作*)
AUOperationRes={{PT,PTk}};(*整合所有的反幺正操作*)
OperationRes={UOperationRes,AUOperationRes};(*将所有对称操作整合起来,总体上分为两类,幺正操作和非幺正操作*)

HamKpGet[orderk_,operation_]:=Module[{nk=orderk,Vk,varlist,matdim,numvar,matcon,hammat,hammatC,diffham,diffhamC,diffcoff,diffcoffC,constrain,constrainAU,constrainAll,constrainC,constrainU,tdc},
nk=orderk;(*哈密顿量多项式部分展开的阶数*)
Vk={kx,ky,kz};
varlist=MonomialList[Expand[(((1+Total[Vk]))^nk)]/._Integer*b_->b,"NegativeDegreeReverseLexicographic"];(*展开多项式并只提取多项式部分，扔掉系数*)
matdim=Length[operation[[1,1,1]]];(*想要展开矩阵的维度,其实就是看表示的维度是多大*)
numvar=Length[varlist];(*哈密顿量展开,多项式的数量*)
tdc=Array[ci,{(Sum[i,{i,matdim}]+Sum[i,{i,matdim-1}])*numvar}];
matcon[numvar_,matdim_]:=Module[{mdr,mdi},mdr=Block[{ta,im},ta=Table[IdentityMatrix[matdim],{j,1,numvar}];(*根据表示的维度构建单位矩阵，数量与多项式的数目相同*)im=0;
Do[Do[im=im+1;
ta[[i,ii,ij]]=tdc[[im]],{ii,1,matdim},{ij,ii,matdim}],{i,1,numvar}];(*i索引ta保存的每一个单位矩阵，ii和ij则是单位矩阵中元素的索引*)Do[Do[ta[[j,ij,ii]]=ta[[j,ii,ij]],{ii,1,matdim},{ij,ii,matdim}],{j,1,numvar}];(*这里是对矩阵中的元素做出一些限制,也就是令ij位置元素和ji位置元素相同，即对称矩阵要求ij=ji*)ta];
mdi=Block[{ta,im},ta=Table[IdentityMatrix[matdim],{j,1,numvar}];(*这里根据多项式展开的数量来确定需要构建的单位矩阵的数量*)im=Sum[i,{i,matdim}]*numvar;(*这里就是构建一个对称矩阵和反对称矩阵来表示厄米的哈密顿量，而自然这两种矩阵中的参数是要不相同的，所以这里就是用im这个索引来区分这些参数，因此在这里的时候使用到的所以就是在上面一个矩阵索引的基础上继续进行叠加，所以这里也就是为什么tdc要构建的dc的数量是(Sum[i,{i,id}]+Sum[i,{i,id-1}])*imd个*)Do[Do[If[ii==ij,ta[[i,ii,ij]]=0,im=im+1;
ta[[i,ii,ij]]=I*tdc[[im]]],{ii,1,matdim},{ij,ii,matdim}],{i,1,numvar}];(*此时限制对角元素为零，非对角和mdr相同只不过乘Parall以了虚数I*)Do[Do[ta[[j,ij,ii]]=-ta[[j,ii,ij]],{ii,1,matdim},{ij,ii,matdim}],{j,1,numvar}];(*此时给出的是限制矩阵为反对称矩阵:ij=-ji*)ta];{mdr,mdi}];
hammat[numvar_,matdim_]:=Module[{mat1},mat1=matcon[numvar,matdim];
varlist.(mat1[[1]]+mat1[[2]])];(*对称矩阵*)
hammatC[numvar_,matdim_]:=Module[{mat1},mat1=matcon[numvar,matdim];
varlist.(mat1[[1]]-mat1[[2]])];(*反对称矩阵构建*)
diffham[opmat_,opk_,numvar_,matdim_]:=opmat.hammat[numvar,matdim].Inverse[opmat]-(hammat[numvar,matdim]/.opk);
diffhamC[opmat_,opk_,numvar_,matdim_]:=opmat.hammatC[numvar,matdim].Inverse[opmat]-(hammat[numvar,matdim]/.opk);
diffcoff[opmat_,opk_,numvar_,matdim_]:=Flatten[CoefficientList[diffham[opmat,opk,numvar,matdim],Vk]];
diffcoffC[opmat_,opk_,numvar_,matdim_]:=Flatten[CoefficientList[diffhamC[opmat,opk,numvar,matdim],Vk]];
constrain[opmat_,opk_,numvar_,matdim_]:=Module[{caplist(*这个参数我是想作为方程求解的时候的变量,因为通过对称性限制之后,一些参数可能会被限制为零,那么在方程求解过程中这些参量就不会作为变量来求解了*),opcons(*这个参量存储对称限制后相差的部分*),re1(*这个变量就用来存储系数方程的解*)},
opcons=diffcoff[opmat,opk,numvar,matdim];
re1=Flatten[Quiet@Solve[opcons==Table[0,{i,1,Length[opcons]}]]];tdc=(tdc/.re1)];


constrainC[opmat_,opk_,numvar_,matdim_]:=Module[{caplist(*这个参数我是想作为方程求解的时候的变量,因为通过对称性限制之后,一些参数可能会被限制为零,那么在方程求解过程中这些参量就不会作为变量来求解了*),opcons(*这个参量存储对称限制后相差的部分*),re1(*这个变量就用来存储系数方程的解*)},
opcons=diffcoffC[opmat,opk,numvar,matdim];
re1=Flatten[Quiet@Solve[opcons==Table[0,{i,1,Length[opcons]}]]];tdc=(tdc/.re1)];

constrainU[OpRes_,numvar_,matdim_]:=Module[{coff},Do[coff=constrain[OpRes[[i,1]],OpRes[[i,2]],numvar,matdim],{i,1,Length[OpRes]}];
coff];(*遍历所有的幺正操作操作进行约束*)


constrainAU[OpRes_,numvar_,matdim_]:=Module[{coff},Do[coff=constrainC[OpRes[[i,1]],OpRes[[i,2]],numvar,matdim],{i,1,Length[OpRes]}];
coff];(*遍历所有的反幺正操作操作进行约束*)


constrainAll[OpRes_,numvar_,matdim_]:=Module[{},constrainU[OpRes[[1]],numvar,matdim];constrainAU[OpRes[[2]],numvar,matdim];];(*整合幺正与非幺正操作*)


constrainAll[OperationRes,numvar,matdim];
tdc=Flatten[tdc/.Quiet[Solve[Variables[tdc]==Table[(Subscript[c, i]),{i,1,Length[Variables[tdc]]}],Variables[tdc]]]];(*这里看起来是个方程求解,其实就是想将之前的构建矩阵时候的参量换个形式表示出来*)
Style[MatrixForm[FullSimplify[hammat[numvar,matdim]/.{kx->(Subscript[k, x]),ky->(Subscript[k, y]),kz->(Subscript[k, z])}]],FontFamily->"Times",20,Blue]
]
HamKpGet[1,OperationRes]