Automatic calculation of generators

zhangzeyingvv · Feb 17, 2023 · 6bee2b6 · 6bee2b6
1 parent 86ec5cf
commit 6bee2b6
Show file tree

Hide file tree

Showing 2 changed files with 45 additions and 8 deletions.
diff --git a/MagneticTB/Utilities.wl b/MagneticTB/Utilities.wl
@@ -143,6 +143,48 @@ Grid[grid,Frame->All]
 ]
 
 
+
+GenerateGroup[gens_,identityElement_,multiply_]:=Module[{i,j,ng,MAXORDER=200,orders,subs,mlist,g1,g2,g3},ng=Length[gens];orders=subs=Range[ng]*0;
+For[i=1,i<=ng,i++,mlist=FoldList[multiply,Table[gens[[i]],MAXORDER]];
+orders[[i]]=FirstPosition[mlist,identityElement][[1]];
+subs[[i]]=mlist[[;;orders[[i]]]];];
+g1=Union@@subs;
+g2=Union@@Table[multiply[g1[[i]],g1[[j]]],{i,Length[g1]},{j,Length[g1]}];
+g3=Complement[g2,g1];g1=g2;
+While[g3!={},g2=Union@@Table[multiply[g3[[i]],g1[[j]]],{i,Length[g3]},{j,Length[g1]}];
+g3=Complement[g2,g1];g1=g2;];
+g2=SortBy[g2,#!=identityElement&]];
+getGenerator[groupele_,identityElement_,multiply_]:=Module[
+{try,group,tmptry,tmpgroup,seteq,groupeleDisorder,gereratorlist},
+(*Greedy Alg, may not give the minimal Generator set*)
+seteq[s1_,s2_]:=SubsetQ[s1,s2]&&SubsetQ[s2,s1];
+If[groupele=={identityElement},Return[groupele]];
+gereratorlist=Table[try={};
+(*groupeleDisorder=RandomSample[groupele];*)
+groupeleDisorder=groupele;
+Do[If[ele==identityElement,Continue[]];
+group=GenerateGroup[try,identityElement,multiply];
+tmptry=Append[try,ele];
+(*Print[group];*)
+tmpgroup=GenerateGroup[tmptry,identityElement,multiply];
+If[Not@seteq[group,tmpgroup],try=tmptry];
+If[seteq[groupele,group],Break[]];,{ele,groupeleDisorder}];
+try,1];
+First[SortBy[gereratorlist,Length[#]&]]];
+findgenind[sgop_]:=Module[{opm,times,identityElement,gens},
+(*opm=MapAt[Mod[#,1]&,sgop[[;;,2;;]],{;;,2}];
+opm=MapAt[#/.{"T"->1,"F"->0}&,opm,{;;,3}];*)
+opm=MapAt[#/.{"T"->1,"F"->0}&,sgop,{;;,4}];
+opm=opm[[;;,{2,4}]];
+(*times=FullSimplify@{#1[[1]].#2[[1]],Mod[#1[[1]].#2[[2]]+#1[[2]],1],Mod[(#1[[3]]+#2[[3]]),2]}&;
+identityElement={IdentityMatrix[3],{0,0,0},0};*)
+times=(*Full*)Simplify@{#1[[1]] . #2[[1]](*,Mod[#1[[1]].#2[[2]]+#1[[2]],1]*),Mod[(#1[[2]]+#2[[2]]),2]}&;
+identityElement={IdentityMatrix[3](*,{0,0,0}*),0};
+gens=getGenerator[opm,identityElement,times];
+Flatten[FirstPosition[opm,#]&/@gens]
+]
+
+
 End[]
 EndPackage[]
 
diff --git a/MagneticTB/magtb.wl b/MagneticTB/magtb.wl
@@ -84,7 +84,8 @@ init[OptionsPattern[]]:=Module[{norm,repall,symopinit},
   nbasesperwyck=Length/@basis;
   Flatten[Table[Table[Table[atompos[[i,j,1]],nbasesperwyck[[i]]],{j,natomperwyck[[i]]}],{i,Length[natomperwyck]}],2]
   ];
-
+  genindex=findgenind[symminfo];
+  Print["Generators:" ,symminfo[[;;,{1,4}]][[genindex]]];
 
 ];
 
@@ -274,12 +275,6 @@ nbases,bases,pbase,pbases,coeff,opmatrix,opmatrixs,tranU,spinrule,solve
     ,{rule,opsrule}]
   ,Length[bases[[1]]]==2,
   Print["Double Value Rep."];
-(*  spinMatrix[op_] := 
-    FullSimplify[Module[{\[Alpha], \[Beta], \[Gamma], nop},
-      If[Det[op] == -1, nop = -op, nop = op];
-        {\[Alpha], \[Beta], \[Gamma]} = -EulerAngles[nop];
-        Transpose@{{E^(I \[Gamma]/2) Cos[\[Beta]/2] E^(I \[Alpha] /2), E^(I \[Gamma]/2) Sin[\[Beta]/2] E^(-I \[Alpha] /2)},
-        {-E^(-I \[Gamma]/2) Sin[\[Beta]/2] E^(I \[Alpha] /2), E^(-I \[Gamma]/2) Cos[\[Beta]/2] E^(-I \[Alpha] /2)}}]];*)
   opsrule=FullSimplify@Table[{MapThread[Rule,{{x,y,z},(Inverse[symm[[2]]] . ({x, y, z} . Inverse[latt]) . latt)}],symm[[4]]},{symm,ops}];
   (*Table[Print[{Transpose[latt] . symm[[2]] . Inverse@Transpose@latt,spinMatrix2[Transpose[latt] . symm[[2]] . Inverse@Transpose@latt]}],{symm,ops}];*)
   spinrule=Table[ExpToTrig@spinMatrix2[Transpose[latt] . symm[[2]] . Inverse@Transpose@latt],{symm,ops}];
@@ -447,7 +442,7 @@ unsymham[n_]:=Module[{bondends,bandind,hop,bonds,toband,nbh,hami},
 ];
 
 
-Options[symham]={symmetryset->All};
+Options[symham]={symmetryset->genindex};
 (*SetOptions[symham,symmetryset\[Rule]All];*)
 symham[n_,OptionsPattern[]]:=Module[{conjh,h,phpmh,recR,para,sset,opset},
   sset=OptionValue[symmetryset];