(*lattice primitive vectors*)
latticeVec={{0,1/2,1/2},{1/2,0,1/2},{1/2,1/2,0}}
(*reciprocal lattice primitive vectors*)
reciprocalVec=Transpose[Inverse[latticeVec]]
(*lattice constant for Germanium*)
aLattice=10.69;
(*half of the volume of the primitive cell*)
vAtomic=(1/8) aLattice^3;

(*important points in the reciprocal space*)
Γ={0.,0.,0.};
X={1.,0.,0.};
L={(1./2.),(1./2.),(1./2.)};
(*basis atoms*)
rb={{-(1./8.),-(1./8.),-(1./8.)},{(1./8.),(1./8.),(1./8.)}};
(*number of atoms in the basis*)
nat=2;

(*local pseudopotential*)
Ra=0.94;
ktf=0.52;
qz=1.83;
Vl[k_]:=If[Abs[k]<(10^(-12)),0,N[16π](((k^2)/(qz^2))-1)((E^(-(((k^2)(Ra^2))/4.)))/((k^2)+(ktf^2)))]

(*non local pseudopotential*)
constnl=12;
Rb=0.71;
lmaxnl=0;
fnl[l_,k_]:=constnl(((Sqrt[π])(Rb^3))/4)If[ l==0,1,(((((Rb k)/2)))^l)](E^(-(((Rb^2) (k^2))/4)))
θ[x_,y_,z_]:=If[(Sqrt[(x^2)+(y^2)+(z^2)])<(10^(-12)),0,ArcCos[(z/(Sqrt[(x^2)+(y^2)+(z^2)]))]];
ϕ[x_,y_]:=Arg[x+I y];
anl[l_,m_,x_,y_,z_,kx_,ky_,kz_]:= 4 π (I^l)(E^(I 2π(x kx+y ky+z kz))) (Sqrt[(((2l+1))/(4 π))])SphericalHarmonicY[l,m,θ[kx,ky,kz],ϕ[kx,ky]]fnl[l,((2π)/aLattice)(Sqrt[(kx^2)+(ky^2)+(kz^2)])];

bso[l_,j_]:=If[j<l,1,-1]

(*definition of vectors in the reciprocal space*)
nsize=5;
indicesTri=Flatten[Outer[List,Table[i,{i,-nsize,nsize}],Table[i,{i,-nsize,nsize}],Table[i,{i,-nsize,nsize}]],2];
indicesDuoTri=Table[{0,0,0},{i,1,2Dimensions[indicesTri][[1]]}];
Do[
indicesDuoTri[[2i-1]]=indicesTri[[i]];
indicesDuoTri[[2i]]=indicesTri[[i]];
,{i,1,Dimensions[indicesTri][[1]]}]
gspaceVec=indicesDuoTri.reciprocalVec;

eMax=10.0; (*maximum kinetic energy*)
kmin=0;
kmax=1;
p=24; (*number of points calculated in one direction*)
nBand=12;(*number of bands to be drawn*)

EnLΓ=Table[0,{x,0,(p/2)}];

(*cycle to calculate the band structure*)
Do[

(*vector k*)
veck=L+2(Γ-L)(x/p);

(*all the k+G*)
kplusgVec=ParallelTable[veck,{i,Dimensions[gspaceVec][[1]]}]+N[gspaceVec];

(*ordering the k+G by increasing energy*)
kplusgVecOrdered=Sort[kplusgVec,(#1.#1<#2.#2)&];

(*select the k+G for which the kinetic energy is<eMax*)
kplusgVecSelected=Select[kplusgVecOrdered,(0.5*(N[2 Pi]/aLattice)^2*#.#<eMax)&];

(*kinetic energy operator*)
kineticOperator=DiagonalMatrix[ParallelTable[(1/2)*(N[2 Pi]/aLattice)^2*(kplusgVecSelected[[i]].kplusgVecSelected[[i]]),{i,Dimensions[kplusgVecSelected][[1]]}]];

(*local pseudopotential operator*)localOperator=Table[If[EvenQ[j+k],(1/(2vAtomic))*Vl[(N[2 Pi]/aLattice)*Sqrt[(kplusgVecSelected[[k,1]]-kplusgVecSelected[[j,1]])^2+(kplusgVecSelected[[k,2]]-kplusgVecSelected[[j,2]])^2+(kplusgVecSelected[[k,3]]-kplusgVecSelected[[j,3]])^2]]*(∑_(n=1)^nat)(Exp[I 2 N[Pi] ((kplusgVecSelected[[k,1]]-kplusgVecSelected[[j,1]])rb[[n]][[1]]+(kplusgVecSelected[[k,2]]-kplusgVecSelected[[j,2]])rb[[n]][[2]]+(kplusgVecSelected[[k,3]]-kplusgVecSelected[[j,3]])rb[[n]][[3]])]),0],{k,Dimensions[kplusgVecSelected][[1]]},{j,Dimensions[kplusgVecSelected][[1]]}];

(*non local projectors*)
anlOperator=Table[anl[l,m,rb[[n]][[1]],rb[[n]][[2]],rb[[n]][[3]],kplusgVecSelected[[i,1]],kplusgVecSelected[[i,2]],kplusgVecSelected[[i,3]]],{n,1,nat},{l,0,lmaxnl},{m,-l,l},{i,Dimensions[kplusgVecSelected][[1]]}];

conjanlOperator=Conjugate[anlOperator];

(*non local operator*)
VKBnl=ParallelTable[If[EvenQ[j+k],(1/(2vAtomic))(∑_(n=1)^natC)(∑_(l=0)^lmaxnl)(∑_(m=-l)^l)(conjanlOperator[[n]][[l+1]][[m+l+1]][[k]]*anlOperator[[n]][[l+1]][[m+l+1]][[j]]),0],{k,Dimensions[kplusgVecSelected][[1]]},{j,Dimensions[kplusgVecSelected][[1]]}];

(*Diagonalize the Hamiltonian and order the eigenvalues by increasing order of energy*)

EnLΓ[[x+1]]=Take[Sort[Eigenvalues[kineticOperator+localOperator+VKBnl]],2nBand];

,{x,0,(p/2)}]

LΓ=Table[{-(p/2)+x,Re[EnLΓ[[x+1]][[j]]-EnLΓ[[p/2+1]][[8]]]*27.211},{j,1,2nBand},{x,0,(p/2)}];

EnΓX=Table[0,{x,0,p}];

(*cycle to calculate the band structure*)
Do[

(*vector k*)
veck=Γ+(X-Γ)(x/p);

(*all the k+G*)
kplusgVec=ParallelTable[veck,{i,Dimensions[gspaceVec][[1]]}]+N[gspaceVec];

(*ordering the k+G by increasing energy*)
kplusgVecOrdered=Sort[kplusgVec,(#1.#1<#2.#2)&];

(*select the k+G for which the kinetic energy is<eMax*)
kplusgVecSelected=Select[kplusgVecOrdered,(0.5*(N[2 Pi]/aLattice)^2*#.#<eMax)&];

(*kinetic energy operator*)
kineticOperator=DiagonalMatrix[ParallelTable[(1/2)*(N[2 Pi]/aLattice)^2*(kplusgVecSelected[[i]].kplusgVecSelected[[i]]),{i,Dimensions[kplusgVecSelected][[1]]}]];

(*local pseudopotential operator*)localOperator=Table[If[EvenQ[j+k],(1/(2vAtomic))*Vl[(N[2 Pi]/aLattice)*Sqrt[(kplusgVecSelected[[k,1]]-kplusgVecSelected[[j,1]])^2+(kplusgVecSelected[[k,2]]-kplusgVecSelected[[j,2]])^2+(kplusgVecSelected[[k,3]]-kplusgVecSelected[[j,3]])^2]]*(∑_(n=1)^nat)(Exp[I 2 N[Pi] ((kplusgVecSelected[[k,1]]-kplusgVecSelected[[j,1]])rb[[n]][[1]]+(kplusgVecSelected[[k,2]]-kplusgVecSelected[[j,2]])rb[[n]][[2]]+(kplusgVecSelected[[k,3]]-kplusgVecSelected[[j,3]])rb[[n]][[3]])]),0],{k,Dimensions[kplusgVecSelected][[1]]},{j,Dimensions[kplusgVecSelected][[1]]}];

(*non local projectors*)
anlOperator=Table[anl[l,m,rb[[n]][[1]],rb[[n]][[2]],rb[[n]][[3]],kplusgVecSelected[[i,1]],kplusgVecSelected[[i,2]],kplusgVecSelected[[i,3]]],{n,1,nat},{l,0,lmaxnl},{m,-l,l},{i,Dimensions[kplusgVecSelected][[1]]}];

conjanlOperator=Conjugate[anlOperator];

(*non local operator*)
VKBnl=ParallelTable[If[EvenQ[j+k],(1/(2vAtomic))(∑_(n=1)^natC)(∑_(l=0)^lmaxnl)(∑_(m=-l)^l)(conjanlOperator[[n]][[l+1]][[m+l+1]][[k]]*anlOperator[[n]][[l+1]][[m+l+1]][[j]]),0],{k,Dimensions[kplusgVecSelected][[1]]},{j,Dimensions[kplusgVecSelected][[1]]}];

(*Diagonalize the Hamiltonian and order the eigenvalues by increasing order of energy*)
EnΓX[[x+1]]=Take[Sort[Eigenvalues[kineticOperator+localOperator+VKBnl]],2nBand];

,{x,0,p}]

ΓX=Table[{x/2,Re[EnΓX[[x+1]][[j]]-EnΓX[[1]][[8]]]*27.211},{j,1,2nBand},{x,0,p}];

(*band structure*)
ListPlot[Flatten[Table[{LΓ[[i]],ΓX[[i]]},{i,1,2nBand}],1], PlotRange->All,Axes->False,Frame->True,FrameLabel->{"k","E"},GridLines->Automatic,GridLinesStyle->Dashed,PlotLegends->Flatten[Table[{"LΓs1"i,"ΓXs1"i,"LΓs2"i,"ΓXs2"i},{i,1,nBand}],1]]

(*Transitions in Γ*)

(*3.0*)
LΓ[[12]][[p/2+1]][[2]]-LΓ[[8]][[p/2+1]][[2]]

(*12.9 PRB32*)
LΓ[[8]][[p/2+1]][[2]]-LΓ[[1]][[p/2+1]][[2]]

(*Transitions in em L*)

(*7.9 PRB33*)
LΓ[[16]][[1]][[2]]-LΓ[[8]][[p/2+1]][[2]]

(*4.2 PRB33*)
LΓ[[12]][[1]][[2]]-LΓ[[8]][[p/2+1]][[2]]

(*0.74 LB*)
LΓ[[10]][[1]][[2]]-LΓ[[8]][[p/2+1]][[2]]

(*1.4*)
LΓ[[8]][[p/2+1]][[2]]-LΓ[[8]][[1]][[2]]

(*7.7*)
LΓ[[8]][[p/2+1]][[2]]-LΓ[[4]][[1]][[2]]

(*10.6*)
LΓ[[8]][[p/2+1]][[2]]-LΓ[[1]][[1]][[2]]