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

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

(*important points in the reciprocal space*)
Γ={0,0,0};
X={1,0,0};
L={(1/2),(1/2),(1/2)};

eMax=6.0;(*maximum kinetic energy*)
p=200; (*number of points calculated in one direction*)

(*cycles to calculate the band structure*)

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

nBand=12; (*number of bands to be drawn*)
Do[
(*vector k*)
kVec=L+2(Γ-L)(x/p);

(*all the k+G*)kplusgVec=Table[kVec,{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,((1/2)*(N[2 Pi]/aLattice)^2*#.#<eMax)&];

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

(*eigenvalues of energy*)
EnLΓ[[x+1]]=Take[Sort[Eigenvalues[kineticOperator]],nBand];

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

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

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

Do[
(*vector k*)
kVec=Γ+(X-Γ)(x/p);

(*all the k+G*)kplusgVec=Table[kVec,{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,((1/2)*(N[2 Pi]/aLattice)^2*#.#<eMax)&];

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

(*eigenvalues of energy*)
EnΓX[[x+1]]=Take[Sort[Eigenvalues[kineticOperator]],nBand];

,{x,0,p}]

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

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