% Tunneling magnetoresistance in 1D model of magnetic tunnel junctions (MTJs)
% Device setup:  F / I / F  <=> ... f f f / o o o o o / f f f ...

clear all; 

%Parameters:
N_b=5; %length of the tunnel barrier region
N_m=2*N_b; %size of the Hamiltonian and Green function matrices
t=1.0; %nearest neighbor hopping in the tunnel barrier region
t_lead=1.0*t; %hooping in the F leads
t_couple=1.0*t; %hopping between the F leads and the tunnel barrier region
E_b=-2*t; %bottom of the band for 1D tight-binding chain 
fermi_energy=0.5*t; %Fermi energy is really measured from the bottom of the band: E_F=fermi_energy-E_b
delta=0.9*(fermi_energy-E_b); %exchange splitting is some fraction of E_F in the ferromagnetic electrodes

%Declare variables explicitly by setting them to zero 
%(not necessary in Matalb but good practice required in C, Fortran90, etc.)
gr=zeros(N_m,N_m); %declare retarded Green function

%Self-energy Sigma of each F electrode is 2x2 matrix in spin space: 
%       (Sigma_up    0      )
%       (   0     Sigma_down)
sigma_left=zeros(2,2); %self-energy due to left F electrode initialized
sigma_right=zeros(2,2); %self-energy due to right F electrode initialized

%Sigma_up: 
energy=fermi_energy-(-delta/2); %exchange splitting for spin-up electrons is -delta/2
rad=4*t_lead^2-energy^2;
if (rad>0.0)
        sigma_up=(energy-i*sqrt(rad));
    else
        sigma_up=(energy-sign(energy)*sqrt(-rad));
end
sigma_up=sigma_up*t_couple^2/(2*t_lead^2);

%Sigma_down:
energy=fermi_energy-(delta/2); %exchange splitting for spin-down electrons is delta/2
rad=4*t_lead^2-energy^2;
if (rad>0.0)
        sigma_down=(energy-i*sqrt(rad));
    else
        sigma_down=(energy-sign(energy)*sqrt(-rad));
end
sigma_down=sigma_down*t_couple^2/(2*t_lead^2);  

point=1; %start counting the number of points in the plot
for U_b=0:0.01:6 %loop over a range of tunnel barrier potential heights (in units of t)
 
barrier(point)=U_b; %vector containing values of the barrier potential

%set up Hamiltonian matrix:
hamil=zeros(N_b,N_b); %without spin, Hamiltonian is N_b x N_b matrix
onsite_potential=U_b-delta/2; %on site potential in the tunnel barrier region where bottom of the band is aligned with spin-up
hamil(1,1)=onsite_potential;
hamil(N_b,N_b)=onsite_potential;
for j=2:N_b %loop over sites
   hamil(j,j)=onsite_potential;
   hamil(j,j-1)=-t;
   hamil(j-1,j)=-t;
 end %loop over sites
 hamil=kron(hamil,eye(2)); %Hamiltonian becomes 2Nb x 2Nb matrix in the Hilbert space H_orbital x H_spin

 left_edge_potential=hamil(1:2,1:2); %record the Hamiltonian matrix element onto which the left self-energy will be added
 right_edge_potential=hamil(N_m-1:N_m,N_m-1:N_m); %record the Hamiltonian matrix element onto which the right self-energy will be added

 %compute transmission function T_P for paralell configuration of magnetizations:
 sigma_left=[sigma_up 0; 0 sigma_down]; %self-energy due to left F electrode
 sigma_right=[sigma_up 0; 0 sigma_down]; %self-energy due to right F electrode
 hamil(1:2,1:2)=left_edge_potential + sigma_left; %add sigma_left to hamiltonian
 hamil(N_m-1:N_m,N_m-1:N_m)=right_edge_potential + sigma_right; %add sigma_right to hamiltonian  
 gr=inv(fermi_energy*eye(N_m)-hamil); %retarded Green function is obtained by brute force inversion of small E-H-sigma_left-sigma_right matrix
 T_parallel=transmission(sigma_left,sigma_right,gr); %call transmission function to use Caroli et al. formula
 T_P(point)=T_parallel;

 %compute transmission function T_AP for anti-paralell configuration of magnetizations:
 sigma_left=[sigma_up 0; 0 sigma_down]; %self-energy due to left F electrode
 sigma_right=[sigma_down 0; 0 sigma_up]; %self-energy due to right F electrode
 hamil(1:2,1:2)=left_edge_potential + sigma_left; %add sigma_left to hamiltonian
 hamil(N_m-1:N_m,N_m-1:N_m)=right_edge_potential + sigma_right; %add sigma_right to hamiltonian
 gr=inv(fermi_energy*eye(N_m)-hamil); %retarded Green function is obtained by brute force inversion of small E-H-sigma_left-sigma_right matrix
 T_antiparallel=transmission(sigma_left,sigma_right,gr); %call transmission function to use Caroli et al. formula
 T_AP(point)=T_antiparallel;
 
point=point+1; %count the number of points in the plot
end %loop over tunnel barrier potential

%compute optimistic tunneling magnetoresistance TMR=(T_P - T_AP)/T_AP at each U_b:
TMR=(T_P-T_AP)./T_AP*100; %multiplied by 100 to get final result in procentage

figure(1); % Select the figure
subplot(2,2,1); % Subdivide the figure into 4 quadrants
plot(barrier,TMR);
grid;
title('TMR vs. Tunneling Barrier Height');
xlabel('U_b/t');
ylabel('Optimistic TMR (%)');
subplot(2,2,3); 
plot(barrier,T_P);
grid;
title('Conductance in P configuration');
xlabel('U_b/t');
ylabel('G(e^2/h)');
subplot(2,2,4); 
plot(barrier,T_AP);
grid;
title('Conductance in AP configuration');
xlabel('U_b/t');
ylabel('G(e^2/h)');
