%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALCULATE THE ENERGY LEVELS OF A 1D PARTICLE-IN-THE-BOX WITH FINITE WALLS %
%                                                                           %
% The Schroeding equation is solved by means of a finite differences scheme %
%               J. Planelles and J.I. Climente, September 2011              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% INPUT DATA

Lmin=5.0;       % shortest QD (nm)
Lmax=50.0;      % longest QD (nm)
Lstep=1.0;      % increment of QD size (nm)
Lb_nm=10.0;     % length of the barrier (nm)
m=0.05;         % electron effective mass (times m0)
Vin_eV=0.00;    % potential energy inside the QD (eV)
Vout_eV=0.25;   % potential energy outside the QD (eV)
h_nm=0.25;      % grid step (nm)
neig=4;         % number of states to be calculated (starting from ground state)

% CALCULATION

for L_nm=Lmin:Lstep:Lmax % open loop of L

% 1) Convert into atomic units
kdist=0.0529177; % kdist nm = 1 bohr radius
kener=27.21138;  % kener eV = 1 hartree

L=L_nm/kdist;
Lb=Lb_nm/kdist;
h=h_nm/kdist;
Vin=Vin_eV/kener;
Vout=Vout_eV/kener;

% 2) Define grid

Lt= 2*Lb + L; % total length
n=round(Lt/h) + 1; % number of grid points

% 3) Calculate potential V(x)

for i=1:n
    x(i)= (i-1)*h; % position at i-th grid point
    
    if (x(i) > Lb-1e-12) && (x(i) < Lb+L+1e-12)
        V(i)=Vin;  % inside the QD
    else
        V(i)=Vout; % outside the QD
    end;

end;
    
% 4) Define coefficients of f(i), f(i+/-1)

for i=1:n
    a(i)=1/(m*h^2)+V(i); % coef of f(i) -depends on i
end;
b=-1/(2*m*h^2);          % coef. of f+ and f-


% 5) Build the Hamiltonian matrix

% define the three diagonals (vectors)
dc=a(2:n-1)';    % main diagonal f(i)
di=ones(n-2,1)*b; % subdiagonal f(i-1)
ds=ones(n-2,1)*b; % superdiagonal f(i+1)

% introduce them into a sparse, tridiagonal matrix
matrix=spdiags([di,dc,ds],[-1,0,1],n-2,n-2);

% 6) Diagonalize the matrix
[P,D]=eigs(matrix,neig,'sa');
eners=diag(D)*kener; % grep energy and convert to eV

% WRITE RESULTS TO DISK

% Write energy to external file
file=fopen('energies.dat','a'); % open file to append results

fprintf(file,'%10.1f',L*kdist);   % write QD size (nm)
for i=1:neig
    fprintf(file,'%10.3f',eners(i)*1000); % write eigenvalues (meV)
endfor;
fprintf(file,'\n');

fclose(file);

end; % close loop of L
