Tight-binding Hamiltonian initialization file

This file should be readable by Maple

# comment lines

desc:="Effective Hamiltonian for pi-system of graphene"; optional description

lattice:="gra": lattice code

P:=Matrix(3,2,[ matrix of all site positions in lattice coordinates as columns
2/3,1/3,
1/3,2/3,
 0 , 0 ]):

dir:=<11,13,0>: lattice direction for use in FGS_dpick command

B:=copy(P): optional convenient custom basis for site positions, e.g. making site coordinates positive integers

encode:=proc(mmt)::string: optional custom labeling of TB elements

RL:=8: radius to which the list L is nearly complete

L:=[ labels for irreducible TB elements, initially can be empty
 [1,1,0,0,0]=e0,  preferably ordered by R
 [1,2,1,0,0]=t1,
 [1,2,2,1,0]=t3b,
...
NULL]:

Ladd:=[ optional extension of L meaningless for TB Hamiltonian but useful for model developmnent
 [1,1,7,8,0]=t16d,
 [1,2,9,5,0]=t17i,
 [1,2,7,8,0]=t17f,
...
NULL]:

The following three tables (K, W, tbv) are needed per demand:

K:=[ k-vectors used in W (table of irreducible representations)
 "G"=<0,0,0>,
"DT"=<k,k,-k>,      preferably, provided k-paths should end at provided k-points, e.g. G-DT-H
 "H"=<1/2,1/2,-1/2>,
...
NULL]:

W:=[ irreducible representations as eigenvectors at high-symmetry k-points and k-directions
 "K1" =[<<0,0,1,exp(2*Pi*I/3),exp(-2*Pi*I/3)>>], list(Matrix) whose columns are eigenvectors
 "G1+"=[<<1,1,0,0,0>|<0,0,1,1,1>>],      here is 2D irreducible block
 "G3+"=[<<0,0,1,-1,0>>, <<0,0,1,1,-2>>], here are two 1D irreducible blocks of a doubly degenerate level
...
NULL]:

tbv:=["    K1     G1     G2     M1     M2  ", parameterization of high-symmetry k-points in consistency with W table
"Ref.[1]  -4.13 -11.67   7.16  -6.48  -2.33",
...
NULL]:

The following variables are optional:

tbe:=[ # model=[e0,t1,t2,t3,t4], minmodel=[e0,t1] eigenvalues at high-symmetry k-points and k-directions in a model fittable by these points
 "K1"=e0-3*t2+6*t4            , # doubly degenerate minmodel is a minimal model free from accidental degenaracy
 "G1"=e0+6*t2+6*t4+(3*t1+3*t3),
 "G2"=e0+6*t2+6*t4-(3*t1+3*t3),
 "DT"=[e0+2*t1+2*t2 +(t1+4*t2+3*t3+4*t4)*cos(k) +2*t4*cos(2*k),  -I*(t1+t3)*sin(k), list(string) ordered as 11,21,22,31,32,33,...
       e0-2*t1+2*t2 -(t1-4*t2+3*t3-4*t4)*cos(k) +2*t4*cos(2*k)],
...
NULL]:

Hke:=Matrix(2):  symbolic matrix H(k)
sbs:=seq(seq(seq(op([s*(n+m)*k1+s*n*k2=s*m*k1-s*n*k3,s*n*k1+s*(n+m)*k2=s*m*k2-s*n*k3]),s=[1,-1]),n=1..5),m=0..5):
h11:=e0+2*add(add(cos(n*k||o),o=1..3)*[t2,t4b,t6b][n],n=1..3)
  +2*add(add(cos(n*(k||o-k||(mod1(o+1,3)))),o=1..3)*[t4,t8][n],n=1..2)
  +2*add(add(add(add(`if`(o1=o2,0,cos(n*k||o1-m*k||o2)),o1=1..3),o2=1..3)*[[],[t6]][n,m],m=1..n-1),n=2..2):
h12:=add(add(exp(I*v),v=[k1,0,-k2]+n*[k1-k3,k2-k1,k3-k2])*[t7,t3,t1,t5c][n+3],n=-2..1)
 +add(add(add(exp(I*v),v=subs(k3=-k1-k2,sbs,[k1,k1,0,0,-k2,-k2]
   +n*[k1-k3,k1-k3,k2-k1,k2-k1,k3-k2,k3-k2]
   +m*[k2-k1,k3-k2,k3-k2,k1-k3,k1-k3,k2-k1]))*[[t5b,0],[t3b,t7c],[t5,t7d]][n+2,m],m=1..2),n=-1..1):
Hke[1,1]:=h11:
Hke[2,2]:=h11:
Hke[1,2]:=h12:
Hke[2,1]:=subs(k1=-k1,k2=-k2,k3=-k3,h12):
Hke:=subs(k3=-k1-k2,Hke):

solver:=proc()  symbolic eigenvalue solver for H(k)
local sf,v,h12R,h12I;
sf:=v->collect(simplify(v),L,sort);
h12R:=sf(evalc(Re(Hke[1,2])));
h12I:=sf(evalc(Im(Hke[1,2])));
<Hke[1,1]-sqrt(h12R^2+h12I^2),Hke[1,1]+sqrt(h12R^2+h12I^2)>
end:

haT:=Matrix(3,shape=symmetric):  symbolic hopping amplitudes tensor
haT[1,1]:=a^2*(t1^2+6*t2^2+10*t3^2)/16:
haT[2,2]:=haT[1,1]:
haT[3,3]:=haT[1,1]:

imT:=table([    table of symbolic inverse mass tensors
"G1"=Matrix(3,(o,q)->`if`(o=q,-a^2*(t1+6*t2+10*t3)/8,0),shape=symmetric),
"H2"=Matrix(3,(o,q)->`if`(o=q, a^2*(t1-6*t2+10*t3)/8,0),shape=symmetric)]):