Drawing Band Structure of Silicene Under Strain Using Tight Binding Methods

[anahita]

I have been drawing band structure silicene under strain using tight binding methods.
How to change K points under strain?
Matlab code is corrected?[
MATLAB code for the intended purpose as follows:

close all
clear all
clc
Na=2; % Number of atoms
Nbnd=8; % number of bands
q=0.1;
w=0.30;
aa=2.25*(1+q);
a=3.82*(1+q);
Csoc=0.3500;
EL=0;
POSx(1)=aa*cosd(30);
POSy(1)=aa*sind(30);
POSz(1)=-0.22;
POSx(2)=2*aa*cosd(30);
POSy(2)=2*aa*sind(30);
POSz(2)=0.22;
a1x=a;
a1y=0;
a2x=a*cosd(60);
a2y=a*sind(60);

sho=0;
for is=[0,-1,1]
  for js=[0,-1,1]
  for ks=1:Na
  sho=sho+1;
  X(sho)=POSx(ks)+is*a1x+js*a2x;
  Y(sho)=POSy(ks)+is*a1y+js*a2y;
  Z(sho)=POSz(ks);
  Ax(sho)=is*a1x+js*a2x;  % Vector for uint cell
  Ay(sho)=is*a1y+js*a2y;
  No(sho)=ks;
  L0=sqrt(((POSx(2)-POSx(1))^2+(POSy(2)-POSy(1))^2+(POSz(2)-POSz(1))^2));
  end
  end
end
figure(1)
plot(X,Y,'*')
Ax=Ax/a;
Ay=Ay/a;
fid = fopen('k.txt', 'r');
KP = fscanf(fid, '%g %g', [3 inf]);
KP = KP';
fclose(fid);

kx=(2*pi/2)*KP(:,1);
ky=(2*pi/2)*KP(:,2);

for ik=1:48
H=H0(Nbnd);
for ie=1:4
  H(ie,ie)=H(ie,ie)+POSz(1)*EL;
end
for ie=5:8
  H(ie,ie)=H(ie,ie)+POSz(2)*EL;
end

HSOC=zeros(2*Nbnd);

for is=1:Na
  for js=1:sho
  dis=sqrt(((X(is)-X(js))^2)+((Y(is)-Y(js))^2)+((Z(is)-Z(js))^2));
  if abs(dis-L0)<0.1 & abs(No(is)-No(js))>0
  l=(X(is)-X(js))/dis;
  m=(Y(is)-Y(js))/dis;
  n=(Z(is)-Z(js))/dis;
  if No(is)==1 & No(js)==2
  h1=hamiltonian1(l,m,n);
  h2=hamiltonian2(l,m,n);
  H(1:4,5:8)=H(1:4,5:8)+h1*exp(2i*((kx(ik)*Ax(js)+ky(ik)*Ay(js))))+h2*exp(2i*((kx(ik)*Ax(js)+ky(ik)*Ay(js))));
  end
  if No(is)==2 & No(js)==1
  h1=hamiltonian1(-l,-m,-n);
  h2=hamiltonian2(-l,-m,-n);
  H(5:8,1:4)=H(5:8,1:4)+h1'*exp(2i*((ky(ik)*Ay(js)+kx(ik)*Ax(js))))+h2'*exp(2i*((kx(ik)*Ax(js)+ky(ik)*Ay(js))));
  end
 
  end
  end
  h=pSOC(Nbnd);
  HSOC((is-1)*8+1:is*8,(is-1)*8+1:is*8)=Csoc*h; 
 
end

for ii=1:Nbnd
  for jj=1:Nbnd
  H2(2*ii-1:2*ii,2*jj-1:2*jj)=H(ii,jj)*eye(2);
  end
endE(ik,1:Nbnd)=sort(real(eig(H)));
Ep(ik,1:2*Nbnd)=sort(real(eig(H2+HSOC)));
pl(ik)=(ik-1)/34.4615;
endfor ip=1:8
  plot(pl,(E(:,ip)),'.b')
  hold all
end
hold all
for ik=1:48
for ib=1:Nbnd
  a=Ep(ik,2*ib-1);
  b=Ep(ik,2*ib);
  E1(ik,ib)=max(a,b);
  E2(ik,ib)=min(a,b);
end
end
E1(:,1:Nbnd)=Ep(:,1:2:2*Nbnd);
E2(:,1:Nbnd)=Ep(:,2:2:2*Nbnd);
figure
for ip=1:8
  plot(pl,real(E1(:,ip)),'.b')
  plot(pl,real(E1(:,ip)),'.b')
  hold all
  plot(pl,real(E2(:,ip)),'.b')
  plot(pl,real(E2(:,ip)),'.b')
end

 

[BigJDubs]

hold allxlabel('K-Points')ylabel('Band Structure')title('Silicene Under strain')function h1=hamiltonian1(l,m,n)h1=[0.5*(1+l^2+m^2+n^2), 0.5*(-l^2+m^2+n^2), 0.5*(-l^2-m^2+n^2), -l*m*(1-n^2); 0.5*(-l^2+m^2+n^2), 0.5*(1+l^2-m^2-n^2), -l*m*(1-n^2), 0.5*(-l^2-m^2+n^2); 0.5*(-l^2-m^2+n^2), -l*m*(1-n^2), 0.5*(1+l^2+m^2-n^2), 0.5*(-l^2+m^2+n^2);-l*m*(1-n^2), 0.5*(-l^2-m^2+n^2), 0.5*(-l^2+m^2+n^2), 0.5*(1+l^2-m^2-n^2)];endfunction h2=hamiltonian2(l,m,n)h2=[0.5*(-l^2+m^2-n^2), -l*m*(1+n^2), 0.5*(-l^2-m^2-n^2), 0.5*(1-l^2+m^2+n^2); -l*m*(1+n^2), 0.5*(1-l^2-m^2+n^2), 0.5*(-l^2+m^2-n^2), 0.5*(-l^2-m^2-n^2); 0.5*(-l^2-m^2-n^2), 0.5*(-l^2+m^2