%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize constants   %
%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
UENERGY_IN_JOULE = 1.602192E-19;
ULENGTH_IN_METER = 1.E-10;
UMASS_IN_KG = 1.66053E-27;
UTIME_IN_SECOND = ULENGTH_IN_METER*sqrt(UMASS_IN_KG/UENERGY_IN_JOULE);
UFREQ_IN_HZ = 1./UTIME_IN_SECOND/2/pi;
UFREQ_IN_THZ = UFREQ_IN_HZ / 1E12;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize Si perfect crystal LDOS   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load dos_si.out;
freq = 1:dos_si(2);
Si = sum(dos_si(4:3+dos_si(2)))
clf;
plot((freq-0.5)*dos_si(3)*UFREQ_IN_THZ,dos_si(4:3+dos_si(2))/Si/dos_si(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;

%%%%%%%%%%%%%%%
% 1st layer   %
%%%%%%%%%%%%%%%
% 200 supercell k
load dos_s5h_80.out;
plot(dos_s5h_80(:,1),filtfilt([1./3. 1./3. 1./3.],[1],dos_s5h_80(:,2)),'r-.');

%%%%%%%%%%%%%%%
% 2nd layer   %
%%%%%%%%%%%%%%%
% 200 supercell k
load dos_s5h_155.out;
plot(dos_s5h_155(:,1),filtfilt([1./3. 1./3. 1./3.],[1],dos_s5h_155(:,2)),'g:');

%%%%%%%%%%%%%%%
% Si near Vc  %
%%%%%%%%%%%%%%%
% 200 supercell k
load dos_s5h_200.out;
plot(dos_s5h_200(:,1),filtfilt([1./3. 1./3. 1./3.],[1],dos_s5h_200(:,2)),'b--');

title ('<001> LDOS of Si in 3 layers near Sigma(5) twist GB compared to that in perfect crystal');
axis([0 35 0 0.2]);
text(7.5,0.185,'1st')
text(32.4,0.095,'2nd')
text(16.5,0.18,'3rd')  
text(28.5,0.135,'perfect')


PRINT;

