(*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)))]

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

(*non local pseudopotential*)
constnl=12;
Rb=0.74;
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)])];

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

Vnl[kx_,ky_,kz_]:=(∑_(n=1)^nat)(∑_(l=0)^lmaxnl)(∑_(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]

constso=.279;
Rc=0.71;
lmaxso=2;
fso[l_,j_,k_]:=If[j>l,(Sqrt[(l/(l+1))]),1]constso (((Sqrt[π])(Rc^3))/4) If[ l==0,1,(((((Rc k)/2)))^l)](E^(-(((Rc^2) (k^2))/4)))
aso[l_,j_,mj_,ms_,x_,y_,z_,kx_,ky_,kz_]:=4π (I^l)(E^(I 2π(x kx+y ky+z kz)))(Sqrt[(((2l+1))/(4 π))])SphericalHarmonicY[l,mj-ms,θ[kx,ky,kz],ϕ[kx,ky]]If[(j>l&&mj==j&&ms>0)||(j>l&&mj==-j&&ms<0),1,If[(j>l&&mj==j&&ms<0)||(j>l&&mj==-j&&ms>0),0,(Sqrt[((2l+1+2(j-l)2ms 2mj)/(4l+2))])]]fso[l,j,((2Pi)/aLattice)(Sqrt[(kx^2)+(ky^2)+(kz^2)])]

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

Plot[{fso[0,(1/2),k]*(Sqrt[27.211]),fso[1,(1/2),k]*(Sqrt[27.211]),fso[1,(3/2),k]*(Sqrt[27.211]),fso[2,(3/2),k]*(Sqrt[27.211]),fso[2,(5/2),k]*(Sqrt[27.211])},{k,0,5},PlotRange->All,PlotLegends->{"l=0, j=(1/2)","l=1, j=(1/2)","l=1, j=(3/2)","l=2, j=(3/2)","l=2, j=(5/2)"},Axes->False,Frame->True,FrameLabel->{"k","fso ((Sqrt[eV]))"},GridLines->Automatic,GridLinesStyle->Dashed]

(*spin orbit pseudopotential*)
Vso[kx_,ky_,kz_]:=(∑_(n=1)^nat)(∑_(l=0)^lmaxso)(∑_(j=Abs[l-(1/2)])^(l+(1/2)))(∑_(mj=-j)^j)Conjugate[aso[l,j,mj,(1/2),rb[[n]][[1]],rb[[n]][[2]],rb[[n]][[3]],kx,ky,kz]]*aso[l,j,mj,(1/2),rb[[n]][[1]],rb[[n]][[2]],rb[[n]][[3]],kx,ky,kz]

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

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

(*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=20; (*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)^nat)(∑_(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]]}];

(*spin orbit projectors*)
asoOperator=ParallelTable[aso[l,j,mj,If[EvenQ[i],-(1/2),(1/2)],rb[[n]][[1]],rb[[n]][[2]],rb[[n]][[3]],kplusgVecSelected[[i,1]],kplusgVecSelected[[i,2]],kplusgVecSelected[[i,3]]],{n,1,nat},{l,0,lmaxso},{j,Abs[l-(1/2)],l+(1/2)},{mj,-j,j},{i,Dimensions[kplusgVecSelected][[1]]}];

conjasoOperator=Conjugate[asoOperator];

(*spin orbit operator*)
VKBso=ParallelTable[(∑_(n=1)^nat)(∑_(l=0)^lmaxso)(∑_(j=Abs[l-(1/2)])^(l+(1/2)))(∑_(mj=-j)^j)(conjasoOperator[[n]][[l+1]][[j-Abs[l-(1/2)]+1]][[mj+j+1]][[k]]*bso[l,j]*asoOperator[[n]][[l+1]][[j-Abs[l-(1/2)]+1]][[mj+j+1]][[i]]),{k,Dimensions[kplusgVecSelected][[1]]},{i, Dimensions[kplusgVecSelected][[1]]}];

(*Diagonalize the Hamiltonian and order the eigenvalues by increasing order of energy*)
EnLΓ[[x+1]]=Take[Sort[Eigenvalues[kineticOperator+localOperator+VKBnl+VKBso]],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)^nat)(∑_(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]]}];

(*spin orbit projectors*)
asoOperator=ParallelTable[aso[l,j,mj,If[EvenQ[i],-(1/2),(1/2)],rb[[n]][[1]],rb[[n]][[2]],rb[[n]][[3]],kplusgVecSelected[[i,1]],kplusgVecSelected[[i,2]],kplusgVecSelected[[i,3]]],{n,1,nat},{l,0,lmaxso},{j,Abs[l-(1/2)],l+(1/2)},{mj,-j,j},{i,Dimensions[kplusgVecSelected][[1]]}];

conjasoOperator=Conjugate[asoOperator];

(*spin orbit operator*)
VKBso=ParallelTable[(∑_(n=1)^nat)(∑_(l=0)^lmaxso)(∑_(j=Abs[l-(1/2)])^(l+(1/2)))(∑_(mj=-j)^j)(conjasoOperator[[n]][[l+1]][[j-Abs[l-(1/2)]+1]][[mj+j+1]][[k]]*bso[l,j]*asoOperator[[n]][[l+1]][[j-Abs[l-(1/2)]+1]][[mj+j+1]][[i]]),{k,Dimensions[kplusgVecSelected][[1]]},{i, 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+VKBso]],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}];

(*gráfico bonito*)
ListPlot[Flatten[Table[{LΓ[[i]],ΓX[[i]]},{i,1,2nBand}],1], PlotRange->{-13,9},Axes->False,Frame->True,FrameLabel->{"k","E (eV)"},GridLines->Automatic,GridLinesStyle->Dashed,PlotLegends->Flatten[Table[{"LΓs1"i,"ΓXs1"i,"LΓs2"i,"ΓXs2"i},{i,1,nBand}],1],AspectRatio->{340,360}]

(*Transtions in Γ*)

(*0.296 PRB12*)
LΓ[[8]][[p/2+1]][[2]]-LΓ[[6]][[p/2+1]][[2]]

(*0.888 (PRB12)*)
LΓ[[10]][[p/2+1]][[2]]-LΓ[[8]][[p/2+1]][[2]]

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

(*0.21*)
LΓ[[16]][[p/2+1]][[2]]-LΓ[[14]][[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]]

(*Other transitions*)

(*1.3 LB*)
banda=Interpolation[Take[ΓX[[10]],-p]];
FindMinimum[banda[x],{x,0.85p/2}][[1]]-LΓ[[8]][[p/2+1]][[2]]

(*3.5 PRB32*)
LΓ[[8]][[p/2+1]][[2]]-ΓX[[6]][[p+1]][[2]]

(*9.5 PRB32*)
LΓ[[8]][[p/2+1]][[2]]-ΓX[[4]][[p+1]][[2]]