%Fig1_12a.m %Essential Electron Transport for Device Physics %band structure for bulk Si, Ge, AlAs, %GaAs, InP, or InAs, along the L - Gamma - X - U,K - Gamma symmetry %directions. The program will prompt the user to enter the material to be %modeled. Method is the sp3s* model as described by %Vogl et. al. in J. Phys. Chem. Solids 44, 365 (1983). The parameters %used are given in the paper. clear clc %plotting parameters + fontsizes FS = 12; %label fontsize FSN = 12; %number fontsize LW = 1; %linewidth % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times'); set(0,'DefaultAxesFontSize', FSN); % Change default text fonts. set(0,'DefaultTextFontname', 'Times'); set(0,'DefaultTextFontSize', FSN); %******************************* % Declare Constants and Arrays * %******************************* % Constants eye = complex(0,1); %Imaginary number, sqrt(-1) d1 = [ 1 1 1]/4; %Nearest neighbor directions d2 = [ 1 -1 -1]/4; d3 = [-1 1 -1]/4; d4 = [-1 -1 1]/4; N = 301; % Number of points calculated for each symmetry direction % Parameters for the sp3s* model % Si Ge AlAs GaAs InP InAs soa = [ 0 0 0.4210 0.4210 0.0670 0.4210]/3; soc = [ 0 0 0.0240 0.1740 0.3920 0.3920]/3; Esa = [-4.2000 -5.8800 -7.5273 -8.3431 -8.5274 -9.5381]; Epa = [ 1.7150 1.6100 0.9833 1.0414 0.8735 0.9099]; Esc = [-4.2000 -5.8800 -1.1627 -2.6569 -1.4826 -2.7219]; Epc = [ 1.7150 1.6100 3.5867 3.6686 4.0465 3.7201]; Esea = [ 6.6850 6.3900 7.4833 8.5914 8.2635 7.4099]; Esec = [ 6.6850 6.3900 6.7267 6.7386 7.0665 6.7401]; Vss = [-8.3000 -6.7800 -6.6642 -6.4513 -5.3614 -5.6052]; Vxx = [ 1.7150 1.6100 1.8780 1.9546 1.8801 1.8398]; Vxy = [ 4.5750 4.9000 4.2919 5.0779 4.2324 4.4693]; Vsapc = [ 5.7292 5.4649 5.1106 4.4800 2.2265 3.0354]; Vpasc = [ 5.7292 5.4649 5.4965 5.7839 5.5825 5.4389]; Vseapc = [ 5.3749 5.2191 4.5216 4.8422 3.4623 3.3744]; Vpasec = [ 5.3749 5.2191 4.9950 4.8077 4.4814 3.9097]; %***************** % Get User Input * %***************** mat = menu('Please choose a material.','Si','Ge','AlAs','GaAs','InP','InAs') clc switch mat case 1 matname = 'Si'; case 2 matname = 'Ge'; case 3 matname = 'AlAs'; case 4 matname = 'GaAs'; case 5 matname = 'InP'; case 6 matname = 'InAs'; end uinput = 0; while uinput == 0 direction = input('Please type a 1 for L-\Gamma-X-U,K-\Gamma symmetry directions, or 2 for an arbitrary direction.') if direction == 1 uinput = 1; elseif direction == 2 uinput = 1; else clc disp('Please choose from the options listed.') end end clc if direction == 1 % l, m, and n vectors for Gamma-L, Gamma-X, X-U, and Gamma-K directions, respectively l = [0.5 1 0 1]; m = [0.5 0 1 1]; n = [0.5 0 1 0]; kmax = [2*pi 2*pi 2*pi 1.5*pi]; kmin = [0 0 1.5*pi 0 ]; else % Input l uinput = 0; while uinput == 0 l = input('Please specify l value'); if isnumeric(l) uinput = 1; else clc disp('Please enter a number') end end clc % Input m uinput = 0; while uinput == 0 m = input('Please specify m value') if isnumeric(m) uinput = 1; else clc disp('Please enter a number') end end clc % Input n uinput = 0; while uinput == 0 n = input('Please specify n value') if isnumeric(n) uinput = 1; else clc disp('Please enter a number') end end clc % Input kmax uinput = 0; while uinput == 0 kmax = input('Please specify the maximum value of k') if isnumeric(kmax) uinput = 1; else clc disp('Please enter a number') end end clc % Input kmin uinput = 0; while uinput == 0 kmin = input('Please specify the minimum value of k') if isnumeric(kmin) uinput = 1; else clc disp('Please enter a number') end end clc end %********************** % Execute Main Script * %********************** for d = 1:length(l) % Symmetry direction loop for i = 1:N % Band calculation loop k = [l(d) m(d) n(d)]*(kmin(d) + (kmax(d)-kmin(d))*(N-i)/(N-1)); % Calculate k vector based on symmetry direction p1 = exp(eye*sum(k.*d1)); p2 = exp(eye*sum(k.*d2)); p3 = exp(eye*sum(k.*d3)); p4 = exp(eye*sum(k.*d4)); g0 = (p1+p2+p3+p4)/4; g1 = (p1+p2-p3-p4)/4; g2 = (p1-p2+p3-p4)/4; g3 = (p1-p2-p3+p4)/4; % Create Hamiltonian % |sa> |sc> |pxa> |pya> |pza> |pxc> |pyc> |pzc> |sea> |sec> h = [Esa(mat)/2 Vss(mat)*g0 0 0 0 Vsapc(mat)*g1 Vsapc(mat)*g2 Vsapc(mat)*g3 0 0 ; %|sa> 0 Esc(mat)/2 -Vpasc(mat)*conj(g1) -Vpasc(mat)*conj(g2) -Vpasc(mat)*conj(g3) 0 0 0 0 0 ; %|sc> 0 0 Epa(mat)/2 0 0 Vxx(mat)*g0 Vxy(mat)*g3 Vxy(mat)*g2 0 -Vpasec(mat)*g1; %|pxa> 0 0 0 Epa(mat)/2 0 Vxy(mat)*g3 Vxx(mat)*g0 Vxy(mat)*g1 0 -Vpasec(mat)*g2; %|pya> 0 0 0 0 Epa(mat)/2 Vxy(mat)*g2 Vxy(mat)*g1 Vxx(mat)*g0 0 -Vpasec(mat)*g3; %|pza> 0 0 0 0 0 Epc(mat)/2 0 0 Vseapc(mat)*g1 0 ; %|pxc> 0 0 0 0 0 0 Epc(mat)/2 0 Vseapc(mat)*g2 0 ; %|pyc> 0 0 0 0 0 0 0 Epc(mat)/2 Vseapc(mat)*g3 0 ; %|pzc> 0 0 0 0 0 0 0 0 Esea(mat)/2 0 ; %|sea> 0 0 0 0 0 0 0 0 0 Esec(mat)/2 ]; %|sec> H = [h+h' zeros(10); zeros(10) h+h']; % Compute spin orbital coupling contribution to Hamiltonian hso = zeros(20); hso(3,4) = -eye*soa(mat); hso(3,15) = soa(mat); hso(4,15) = -eye*soa(mat); hso(5,13) = -soa(mat); hso(5,14) = eye*soa(mat); hso(6,7) = -eye*soc(mat); hso(6,18) = soc(mat); hso(7,18) = -eye*soc(mat); hso(8,16) = -soc(mat); hso(8,17) = eye*soc(mat); hso(13,14) = eye*soa(mat); hso(16,17) = eye*soc(mat); Hso = hso + hso'; % Compute eigenvalues of spin orbital modified Hamiltonian [V,D] = eig(H + Hso); eiglst = sum(D); E(i,:,d) = sort(real(eiglst)); end end %*************** % Plot Results * %*************** if direction == 1 % Reorder symmetry directions for plotting and create k-axis % L - Gamma symmetry for i = 1:N plotE(i,:) = E(i,:,1); plotK(i) = -0.5*(N-i)/(N-1); end % Gamma - X symmetry for i = 1:N plotE(N+i,:) = E(N+1-i,:,2); plotK(N+i) = (i-1)/(N-1); end % X - U symmetry for i = 1:N plotE(2*N+i,:) = E(i,:,3); plotK(2*N+i) = 1 + 0.25*(i-1)/(N-1); end % K - Gamma symmetry for i = 1:N plotE(3*N+i,:) = E(i,:,4); plotK(3*N+i) = 1.25 + .75*(i-1)/(N-1); end figure(1) clf plot(plotK,plotE,'LineWidth',LW) title([matname,' dispersion in L-\Gamma-X-U,K-\Gamma symmetry directions']) xlabel({'';'';'';'Normalized Wavenumber, kL'}) ylabel('Energy (eV)') xlim([-0.5 2]) grid on set(gca,'XTickLabel',{''}) text(-0.5,-16,'\pi') text(-0.5,-18,'L') text(0,-16,'0') text(0,-18,'\Gamma') text(1,-16,'2\pi') text(1,-18,'X') text(1.25,-16,'3\pi/2') text(1.25,-18,'U,K') text(2,-16,'0') text(2,-18,'\Gamma') set(gca,'XTick',[-1/sqrt(3) 0 1 1.25 2]) else dK = (kmax-kmin)/(N-1); plotK = kmin:dK:kmax; for i = 1:N plotE(i,:) = E(i,:,1); end plot(plotK,plotE) title(['Energy band diagram for ',matname,' in the [l,m,n]=[',num2str(l),',',num2str(m),',',num2str(n),'] direction with k_m_a_x=',num2str(kmax),' and k_m_i_n=',num2str(kmin),'.']) xlabel('Wavenumber,kL/(2\pi)') ylabel('Energy (eV)') xlim([kmin kmax]) grid on end