Kronig-Penney model matlab code

Back to spectroscopic techniques animations

Download

matlab code

% Program to plot the dispersion curve and density of states of a 1-D array

% of atoms according to the Kronig-Penney model, a simplified model in

% solid-state physics for electrons in a one-dimensional pulselike

% periodic potential. See e.g.,

% https://en.wikipedia.org/wiki/Particle_in_a_one-dimensional_lattice for

% the physics background.

close all

clear

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]);

% Locations and sizes of subplots

pos1 = [0.05 0.25 0.28 0.4];

pos2 = [0.38 0.25 0.28 0.4];

pos3 = [0.71 0.25 0.28 0.4];

nFrames = 300; % Number of frames in the movie

potLength = 8; % Length of potential plot in x-direction in Angstroms

% Constants

hbar = 1.05457E-34; % Reduced Planck's constant

ec = 1.6022E-19; % Elemental charge

% Functions for hyperbolic cosine and sine

cosh = @(x) (exp(x) + exp(-x)) / 2;

sinh = @(x) (exp(x) - exp(-x)) / 2;

prompt = ['Which parameter do you want to scan? \n' ...

'V = potential difference between maximum and minimum in periodic perturbation \n' ...

'P = Periodicity \n' ...

'D = duty cycle (peak width/periodicity): '];

varPar = input(prompt, 's');

prompt = 'Lower energy limit for plot [eV]? Default = -1 eV: ';

E1 = input(prompt);

if isempty(E1)

E1 = -1;

end

prompt = 'Upper energy limit for plot [eV]? Default = 50 eV: ';

E2 = input(prompt);

if isempty(E2)

E2 = 50;

end

prompt = 'Lower potential of periodic square wave [V]? Default = 0 V: ';

V1 = input(prompt);

if isempty(V1)

V1 = 0; % Can't think of a situation where one would not set this to 0

end

% Variable value of pulse potential height

if (varPar == 'v' || varPar == 'V') % Potential is the variable

prompt = 'Starting upper potential of periodic square wave [V]? : ';

V21 = input(prompt);

prompt = 'Ending upper potential of periodic square wave [V]? : ';

V22 = input(prompt);

Vmax = V22;

loopVals = linspace(V21,V22,nFrames);

vid = VideoWriter('kronigPenney_potential.mp4','MPEG-4');

else % Not the variable

prompt = 'Upper potential of periodic square wave [V]? Default = 10 V: ';

V2 = input(prompt);

Vmax = V2;

if isempty(V2)

V2 = 10;

Vmax = V2;

end

% Variable value of pulse potential periodicity

if (varPar == 'p' || varPar == 'P') % Periodicity is the variable

prompt = 'Starting periodicity [Å]? : ';

P1 = input(prompt)*1e-10;

prompt = 'Ending periodicity [Å]? : ';

P2 = input(prompt)*1e-10;

loopVals = linspace(P1,P2,nFrames);

vid = VideoWriter('kronigPenney_periodicity.mp4','MPEG-4');

else % Not the variable

prompt = 'Periodicity of 1-D potential [Å]? Default = 1.6 Å: ';

a = input(prompt);

if isempty(a)

a = 1.6e-10;

else

a = a*1e-10; % Convert to m

end

% Variable value of pulse potential duty cycle

if (varPar == 'd' || varPar == 'D') % Duty cycle is the variable

prompt = 'Starting pulse width [Å]? : ';

c1 = input(prompt)*1e-10;

prompt = 'Ending pulse width [Å]? : ';

c2 = input(prompt)*1e-10;

b1 = a - c1; % Starting width of lower potential in m

b2 = a - c2; % Ending width of lower potential in m

loopVals = linspace(b1,b2,nFrames);

vid = VideoWriter('kronigPenney_dutyCycle.mp4','MPEG-4');

else % Not the variable

prompt = 'Width of peak voltage of 1-D potential [Å]? Default = 0.6 Å: ';

c = input(prompt);

if isempty(c)

c = 0.6e-10;

else

c = c*1e-10; % Convert to m

end

m = 9.11e-31; % Mass of electron in kg

N = 5000; % Number of energy points to plot

vid.Quality = 100;

vid.FrameRate = 30;

open(vid);

% Calculate dispersion and DoS curves

for ii = loopVals

% Initialize arrays

dispersion = zeros(N, 2);

density_of_states = zeros(N, 2);

if (varPar == 'd' || varPar == 'D') % Duty cycle is the variable

b = ii; % Change b-value for each loop

V2 = Vmax;

end

if (varPar == 'v' || varPar == 'V') % Peak potential is the variable

V2 = ii; % Change upper-voltage-value for each loop

b = a - c;

end

if (varPar == 'p' || varPar == 'P') % Periodicity is the variable

a = ii; % Change periodicity for each loop

b = ii - c; % Width of lower potential in m

V2 = Vmax;

end

potLoop = ceil(potLength*1e-10/a);

potential = zeros(4*potLoop,2);

for jj = 0:potLoop-1

potential(1+jj*4,1) = 0 + a*jj*1e9;

potential(1+jj*4,2) = V1;

potential(2+jj*4,1) = b*1e9 + a*jj*1e9;

potential(2+jj*4,2) = V1;

potential(3+jj*4,1) = b*1e9 + a*jj*1e9;

potential(3+jj*4,2) = V2;

potential(4+jj*4,1) = a*1e9 + a*jj*1e9;

potential(4+jj*4,2) = V2;

end

% Convert potentials to Joules

V1 = V1 * ec;

V2 = V2 * ec;

maxDOE = 0.0;

% Loop through energy levels

for i = 1:N % i-loop

E = E1 + (E2 - E1) * (i - 1) / (N - 1);

EJ = E * ec; % energy in Joules

k1 = sqrt(2 * m * abs(EJ - V1)) / hbar;

k2 = sqrt(2 * m * abs(EJ - V2)) / hbar;

dk1dE = m / (hbar^2 * k1);

dk2dE = m / (hbar^2 * k2);

if EJ < V1

alpha = 2 * cosh(k2 * (a - b)) * cosh(k1 * b) + ...

(k2 / k1 + k1 / k2) * sinh(k2 * (a - b)) * sinh(k1 * b);

dispersion(i, :) = [0, E];

density_of_states(i, :) = [0, E];

if (max(density_of_states(i,2)) > maxDOE)

maxDOE = max(density_of_states(i,2));

end

else

if EJ < V2

dk2dE = -m / (hbar^2 * k2);

alpha = 2 * cosh(k2 * (a - b)) * cos(k1 * b) + ...

(k2 / k1 - k1 / k2) * sinh(k2 * (a - b)) * sin(k1 * b);

dalphadk1 = -2 * b * cosh(k2 * (a - b)) * sin(k1 * b) + ...

(-k2 / (k1^2) - 1 / k2) * sinh(k2 * (a - b)) * sin(k1 * b) + ...

b * (k2 / k1 - k1 / k2) * sinh(k2 * (a - b)) * cos(k1 * b);

dalphadk2 = 2 * (a - b) * sinh(k2 * (a - b)) * cos(k1 * b) + ...

(1 / k1 + k1 / (k2^2)) * sinh(k2 * (a - b)) * sin(k1 * b) + ...

(a - b) * (k2 / k1 - k1 / k2) * cosh(k2 * (a - b)) * sin(k1 * b);

else

alpha = 2 * cos(k2 * (a - b)) * cos(k1 * b) - ...

(k2 / k1 + k1 / k2) * sin(k2 * (a - b)) * sin(k1 * b);

dalphadk1 = -2 * b * cos(k2 * (a - b)) * sin(k1 * b) - ...

(-k2 / (k1^2) + 1 / k2) * sin(k2 * (a - b)) * sin(k1 * b) - ...

b * (k2 / k1 + k1 / k2) * sin(k2 * (a - b)) * cos(k1 * b);

dalphadk2 = -2 * (a - b) * sin(k2 * (a - b)) * cos(k1 * b) - ...

(1 / k1 - k1 / (k2^2)) * sin(k2 * (a - b)) * sin(k1 * b) - ...

(a - b) * (k2 / k1 + k1 / k2) * cos(k2 * (a - b)) * sin(k1 * b);

end

dkda = 1 / (a * sqrt(4 - alpha^2));

if alpha > 0

if alpha^2 < 4

dispersion(i, :) = [atan(sqrt(4 - alpha^2) / alpha), E];

density_of_states(i, :) = [abs(ec * 1E-9 * (2 / pi) * ...

dkda * (dalphadk1 * dk1dE + dalphadk2 * dk2dE)), E];

else

density_of_states(i, :) = [0, E];

dispersion(i, :) = [0, E];

end

else

if alpha^2 < 4

dispersion(i, :) = [pi + atan(sqrt(4 - alpha^2) / alpha), E];

density_of_states(i, :) = [abs(ec * 1E-9 * (2 / pi) * ...

dkda * (dalphadk1 * dk1dE + dalphadk2 * dk2dE)), E];

else

density_of_states(i, :) = [0, E];

dispersion(i, :) = [pi, E];

end

if (max(density_of_states(i,2)) > maxDOE)

maxDOE = max(density_of_states(i,2));

end

end % end of i-loop

% Plotting

hold off

subplot('position',pos1); % Plot the rectangular periodic potential

newplot

plot(potential(:, 1), potential(:, 2), 'r','LineWidth', 2.0);

xlabel('Position [nm]');

ylabel('Potential [eV]');

xlim([0 potLength/10]);

ylim([V1,Vmax*1.05]);

set(gca,'fontsize',18);

set(gca,'TickLength',[0.02, 1]);

set(gca,'linewidth',2);

box('off')

if (varPar == 'd' || varPar == 'D') % Duty cycle is the variable

str1 = 'Potential = ';

str2 = num2str(V2/ec,4);

str3 = ' eV';

strTot = [str1 str2 str3];

text(0.0, 1.1, strTot, 'FontSize',18,'units','normalized');

str1 = 'Periodicity = ';

str2 = num2str(a/1e-10,4);

str3 = ' Å';

strTot = [str1 str2 str3];

text(0.35, 1.1, strTot, 'FontSize',18,'units','normalized');

str1 = 'Duty cycle = ';

str2 = num2str((a-ii)/a,4);

strTot = [str1 str2];

text(0.7, 1.1, strTot, 'FontSize',18,'units','normalized');

end

if (varPar == 'v' || varPar == 'V') % Peak potential is the variable

str1 = 'Potential = ';

str2 = num2str(ii,4);

str3 = ' eV';

strTot = [str1 str2 str3];

text(0.0, 1.1, strTot, 'FontSize',18,'units','normalized');

str1 = 'Periodicity = ';

str2 = num2str(a/1e-10,4);

str3 = ' Å';

strTot = [str1 str2 str3];

text(0.37, 1.1, strTot, 'FontSize',18,'units','normalized');

str1 = 'Duty cycle = ';

str2 = num2str((a-b)/a,4);

strTot = [str1 str2];

text(0.7, 1.1, strTot, 'FontSize',18,'units','normalized');

end

if (varPar == 'p' || varPar == 'P') % Periodicity is the variable

str1 = 'Potential = ';

str2 = num2str(V2/ec,4);

str3 = ' eV';

strTot = [str1 str2 str3];

text(0.0, 1.1, strTot, 'FontSize',18,'units','normalized');

str1 = 'Periodicity = ';

str2 = num2str(ii/1e-10,4);

str3 = ' Å';

strTot = [str1 str2 str3];

text(0.32, 1.1, strTot, 'FontSize',18,'units','normalized');

str1 = 'Duty cycle = ';

str2 = num2str((a-b)/a,4);

strTot = [str1 str2];

text(0.7, 1.1, strTot, 'FontSize',18,'units','normalized');

end

hold off

subplot('position',pos2); % Plot the dispersion as a function of ka

newplot

plot(1.001*dispersion(:, 1)-0.001, dispersion(:, 2), 'g','LineWidth', 2.0);

xlabel('ka');

ylabel('E [eV]');

xticks([0 pi/4 pi/2 3*pi/4 pi])

xticklabels({'0','\pi/4','\pi/2','3\pi/4','\pi'})

xlim([0,pi]);

ylim([E1,E2]);

set(gca,'fontsize',18);

set(gca,'TickLength',[0.02, 1]);

set(gca,'linewidth',2);

hold off

subplot('position',pos3); % Plot the density of states as a function of energy

newplot

lastPt = [0 50];

density_of_states = [density_of_states; lastPt];

plot(density_of_states(:, 1), density_of_states(:, 2), 'b','LineWidth', 2.0);

hold on

fill(density_of_states(:, 1), density_of_states(:, 2), 'b','FaceAlpha',0.25);

xlabel('Density of States [eV^{-1}nm^{-1}]');

ylabel('E [eV]');

xlim([0,2.5]);

ylim([E1,E2]);

set(gca,'fontsize',18);

set(gca,'TickLength',[0.02, 1]);

set(gca,'linewidth',2);

frame = getframe(gcf);

writeVideo(vid,frame);

end % loopVals loop

close(vid)