% Fig_1_17.m
% Essential Semiconductor Laser Device Physics
% This program calculates the bulk band structure for Si, Ge, AlAs,
% GaAs, InP, or InAs, along the L - Gamma - X - U,K - Gamma symmetry
% directions.  The program prompts the user to enter the material to be
% modeled.  The method used is the sp3s* model as described by
% Vogl et. al. in J. Phys. Chem. Solids 44, 365 (1983).  The parameters
% used are as given in the paper.

clear all
clf;		
%plotting parameters + fontsizes
FS      = 16;	%label fontsize
FSN     = 14;	%number fontsize
LW      = 2;	%linewidth

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times New Roman');
set(0,'DefaultAxesFontSize', FSN);

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times New Roman');
set(0,'DefaultTextFontSize', FSN);

hbarChar=['\fontname{MT Extra}h\fontname{Times New Roman}'];

%*******************************
% 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('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('Type 1 for L-\Gamma-X-U,K-\Gamma symmetry directions, or 2 for 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('Specify 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('Specify 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
        % Calculate k vector based on symmetry direction
        k = [l(d) m(d) n(d)]*(kmin(d) + (kmax(d)-kmin(d))*(N-i)/(N-1)); 
        
        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 *
%***************
    GLpt=-sqrt(3)/2;%distance from gamma to L point normalized to gamma-X
    GXpt=1;
    XUpt=0.5/sqrt(2);
    GKpt=0.5*3/sqrt(2);
    
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) = GLpt*(N-i)/(N-1);
    end
    
    % Gamma - X symmetry
    for i = 1:N
        plotE(N+i,:) = E(N+1-i,:,2);
        plotK(N+i) = GXpt*(i-1)/(N-1);
    end
    
    % X - U symmetry
    for i = 1:N
        plotE(2*N+i,:) = E(i,:,3);
        plotK(2*N+i) = GXpt + XUpt*(i-1)/(N-1);
    end
    
    % K - Gamma symmetry
    for i = 1:N
        plotE(3*N+i,:) = E(i,:,4);
        plotK(3*N+i) = (GXpt + XUpt) + GKpt*(i-1)/(N-1);
    end

    p0oset=-20;
    p1oset=-18;%offset on x-axis labels
    p2oset=-16;
%********************  plot figure ************************    
    figure(1)
    clf
    plot(plotK,plotE,'LineWidth',LW)
    title(['Energy band diagram for ',matname,...
        ' along the L-\Gamma-X-U,K-\Gamma symmetry directions'])
    %text(0,p0oset,'Normalized Wavenumber, kL');%xlabel
    xlabel({'';'';'';'';'Normalized wave vector, \itkL\rm'})
    ylabel('Energy, \itE\rm (eV)')
    xlim([GLpt (GXpt+XUpt+GKpt)])
    grid on
    set(gca,'XTickLabel',{''})
    text(GLpt*0.6,p1oset,'[111]');
    text(GLpt,p2oset,'L');text(GLpt/2,p2oset,'\Lambda');
    text(0,p2oset,'\Gamma')
    text(GXpt*0.4,p1oset,'[100]');
    text(GXpt,p2oset,'X');text(GXpt/2,p2oset,'\Delta');
    text((GXpt+XUpt),p2oset,'U,K')
    text((GXpt+XUpt+GKpt*0.41),p1oset,'[011]');
    text((GXpt+XUpt+GKpt/2),p2oset,'\Sigma');
    text((GXpt+XUpt+GKpt),p2oset,'\Gamma')
    set(gca,'XTick',[GLpt 0 GXpt (GXpt+XUpt) (GXpt+XUpt+GKpt)])
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('Wave vector, \itkL\rm/(2\pi)')
    ylabel('Energy, \itE\rm (eV)')
    xlim([kmin kmax])
    grid on
end