(*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;

(*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.972;
ktf=0.62;
qz=2.17;
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)))]

Plot[Vl[k]*27.211,{k,0,1},Axes->False,Frame->True,FrameLabel->{"k","Vlocal (eV)"},GridLines->Automatic,GridLinesStyle->Dashed]

(*non local pseudopotential*)
Rb=1.06;
constnl=6.1;
lmax=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) (Sqrt[(((2l+1))/(4 π))])(E^(I 2π(x kx+y ky+z kz)))SphericalHarmonicY[l,m,θ[kx,ky,kz],ϕ[kx,ky]]fnl[l,((2π)/aLattice)(Sqrt[(kx^2)+(ky^2)+(kz^2)])];

Plot[fnl[0,k]*(Sqrt[27.211]),{k,0,6},Axes->False,Frame->True,FrameLabel->{"k","fl ((Sqrt[eV]))"},GridLines->Automatic,GridLinesStyle->Dashed]

bnl[l_]:=Sign[NIntegrate[(x^2)fnl[l,x],{x,0,20}]]
bnl[0]

Vnl[kx_,ky_,kz_]:=(∑_(n=1)^nat)(∑_(l=0)^lmax)(∑_(m=-l)^l)(Conjugate[anl[l,m,rb[[n]][[1]],rb[[n]][[2]],rb[[n]][[3]],kx,ky,kz]]*anl[l,m,rb[[n]][[1]],rb[[n]][[2]],rb[[n]][[3]],kx,ky,kz])

Plot[Vnl[kx,0,0]*27.211,{kx,0,1},PlotRange->All,Axes->False,Frame->True,FrameLabel->{"k","Vnonlocal (eV)"},GridLines->Automatic,GridLinesStyle->Dashed]

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

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

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

Dynamic[x]

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

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

(*all the k+G*)
kplusgVec=Table[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[Table[0.5*(N[2 Pi]/aLattice)^2*(kplusgVecSelected[[i]].kplusgVecSelected[[i]]),{i,Dimensions[kplusgVecSelected][[1]]}]];

(*local pseudopotential operator*)
localOperator=Table[(1./(2.vAtomic))*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]])]),{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,lmax},{m,-l,l},{i,Dimensions[kplusgVecSelected][[1]]}];

conjanlOperator=Conjugate[anlOperator];

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

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

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

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

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

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

(*ciclo para calcular a estrutura de bandas*)
Do[

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

(*all the k+G*)
kplusgVec=Table[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[Table[0.5*(N[2 Pi]/aLattice)^2*(kplusgVecSelected[[i]].kplusgVecSelected[[i]]),{i,Dimensions[kplusgVecSelected][[1]]}]];

(*local pseudopotential operator*)
localOperator=Table[(1./(2.vAtomic))*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]])]),{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,lmax},{m,-l,l},{i,Dimensions[kplusgVecSelected][[1]]}];

conjanlOperator=Conjugate[anlOperator];

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

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

,{x,0,p}]

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

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

(*Transitions of Silicon*)

1.2;
LΓ[[4]][[1]][[2]]

2.04;
LΓ[[5]][[1]][[2]]

3.35;
LΓ[[7]][[p/2+1]][[2]]

4.15;
LΓ[[8]][[p/2+1]][[2]]

banda=Interpolation[ΓX[[5]]]
m=FindMinimum[banda[x],{x,0.85p/2}];

1.17;
m[[1]]

0.85;

(x/.m[[2]])/(p/2)

2.9;
ΓX[[4]][[p+1]][[2]]

12.5;
ΓX[[1]][[1]][[2]]