% 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 ...
function Mtj
clear all;
tic
%Parameters:
d=3;%dimetion of the system. It can be d=1 for 1D chain, d=2 for 2D plane and d=3 for 3D bulk
V_bias=0;% biased voltage for charge and spin claculations
Tmpr=30;%temperature of the system for chrge and spin current calculations
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
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=2.25; %exchange splitting is some fraction of E_F in the ferromagnetic electrodes
de=0.015; %energy grid spacing within the interval -2*E_b to 2*E_b for in-plane 2D tight-binding dispersion
Tmpr=Tmpr/(11604.5);
%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); %#ok<NASGU> %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); %#ok<NASGU> %self-energy due to left F electrode initialized
sigma_right=zeros(2,2); %#ok<NASGU> %self-energy due to right F electrode initialized

%Sigma_up:
 %exchange splitting for spin-up electrons is -delta/2
sigma_up=@(E)self_energy(E-(-delta/2),t_lead);
%Sigma_down:
%exchange splitting for spin-down electrons is delta/2
sigma_down=@(E)self_energy(E-(delta/2),t_lead);

point=1; %start counting the number of points in the plot
U_b=4;
h = waitbar(0,'Initializing...');
for U_b=0:0.02:6 %loop over a range of tunnel barrier potential heights (in units of t)
    waitbar(U_b/6,h,'Progress...');
    barrier(point)=U_b; %vector containing values of the barrier potential
    hamil=Hamiltonian(N_b,U_b-delta/2,t);    
    
    %compute transmission function T_P for paralell configuration of magnetizations:
    sigma_left=@(E)[sigma_up(E) 0; 0 sigma_down(E)]; %self-energy due to left F electrode
    sigma_right=@(E)[sigma_up(E) 0; 0 sigma_down(E)]; %self-energy due to right F electrode
    [I_charge_P(point),I_Sx_P(point),I_Sy_P(point),I_Sz_P(point),T_P(point)]=observables(d,V_bias,t,Tmpr,fermi_energy,hamil,sigma_left,sigma_right);
    
    %compute transmission function T_AP for anti-paralell configuration of magnetizations:
    sigma_left=@(E)[sigma_up(E) 0; 0 sigma_down(E)]; %self-energy due to left F electrode
    sigma_right=@(E)[sigma_down(E) 0; 0 sigma_up(E)]; %self-energy due to right F electrode
    [I_charge_AP(point),I_Sx_AP(point),I_Sy_AP(point),I_Sz_AP(point),T_AP(point)]=observables(d,V_bias,t,Tmpr,fermi_energy,hamil,sigma_left,sigma_right);
    %compute optimistic tunneling magnetoresistance TMR=(T_P - T_AP)/T_AP:
    
    point=point+1; %count the number of points in the plot
end %loop over tunnel barrier potential

close(h)
%Calculate zero bias optimistic TMR
TMR=(T_P-T_AP)./T_AP*100; %multiplied by 100 to get result in procentage
%Calculate finite bias optimistic TMR
TMR_V=(I_charge_P-I_charge_AP)./I_charge_AP*100; %multiplied by 100 to get result in procentage
figure(1); % Select the figure
subplot(2,2,1); % Subdivide the figure into 4 quadrants
plot(barrier,TMR,barrier,TMR_V);
grid;
title('TMR vs. Tunneling Barrier Height');
xlabel('U_b/t');
ylabel('Optimistic TMR (%)');

%plot both transmission and current/V_bias for comparison
subplot(2,2,3);
plot(barrier,T_P,barrier,I_charge_P/V_bias);
grid;
title('Conductance in P configuration');
xlabel('U_b/t');
ylabel('G(e^2/h)');

%plot both transmission and current/V_bias for comparison
subplot(2,2,4);
plot(barrier,T_AP,barrier,I_charge_AP/V_bias);
grid;
title('Conductance in AP configuration');
xlabel('U_b/t');
ylabel('G(e^2/h)');

figure(2); % Select the figure
subplot(1,2,1);
plot(barrier,I_Sx_P/V_bias,barrier,I_Sy_P/V_bias,barrier,I_Sz_P/V_bias);
grid;
title('Spin currents in P configuration');
xlabel('U_b/t');
ylabel('I_S(h/2\pi)');
legend('I_{Sx}','I_{Sy}','I_{Sz}');
subplot(1,2,2);
plot(barrier,I_Sx_AP/V_bias,barrier,I_Sy_AP/V_bias,barrier,I_Sz_AP/V_bias);
grid;
title('Spin currents in AP configuration');
xlabel('U_b/t');
ylabel('I_S(h/2\pi)');
legend('I_{Sx}','I_{Sy}','I_{Sz}');
toc;
end
function SE=self_energy(energy,t)
rad=4*t^2-energy^2;
if (rad>0.0)
    SE=(energy-1i*sqrt(rad));
else
    SE=(energy-sign(energy)*sqrt(-rad));
end
SE=SE*t^2/(2*t^2);
end
function H=Hamiltonian(N,E,t)
    %set up Hamiltonian matrix:
    H=zeros(N,N); %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
    for j=1:N %loop over sites
        H(j,j)=E;
    end %loop over sites
    for j=2:N
        H(j,j-1)=-t;
        H(j-1,j)=-t;
    end
    H=kron(H,eye(2)); %Hamiltonian becomes 2Nb x 2Nb matrix in the Hilbert space H_orbital x H_spin
end
function T=Trans(EF,H,sigma_left,sigma_right)
T=zeros(size(EF));
for i=1:max(size(EF))
    hamil=H;
    N_m=size(hamil,1);
    hamil(1:2,1:2)=hamil(1:2,1:2) + sigma_left(EF(i)); %add sigma_left to hamiltonian
    hamil(N_m-1:N_m,N_m-1:N_m)=hamil(N_m-1:N_m,N_m-1:N_m) + sigma_right(EF(i)); %add sigma_right to hamiltonian
    gr=inv(EF(i)*eye(N_m)-hamil); %retarded Green function is obtained by brute force inversion of small E-H-sigma_left-sigma_right matrix
    T(i)=transmission(sigma_left(EF(i)),sigma_right(EF(i)),gr); %call transmission function to use Caroli et al. formula
end
end
function I=Current(S,mu,Tmpr,Energy,e,H,sigma_left,sigma_right,m)
I=zeros(size(e));
for i=1:max(size(e))
    hamil=H;
    N_m=size(hamil,1);
    Sig1=sigma_left(Energy-e(i));Sig2=sigma_right(Energy-e(i));
    gamma_left=1i*(Sig1-Sig1'); %level broadening matrix from left F electrode
    gamma_right=1i*(Sig2-Sig2'); %level broadening matrix from right F electrode
    F_left=1/(1+exp((Energy-mu(1))/Tmpr));
    F_right=1/(1+exp((Energy-mu(2))/Tmpr));
    
    hamil(1:2,1:2)=hamil(1:2,1:2) + Sig1; %add sigma_left to hamiltonian
    hamil(N_m-1:N_m,N_m-1:N_m)=hamil(N_m-1:N_m,N_m-1:N_m) + Sig2; %add sigma_right to hamiltonian
    gr=inv((Energy-e(i))*eye(N_m)-hamil); %retarded Green function is obtained by brute force inversion of small E-H-sigma_left-sigma_right matrix
    j1=1:2;j2=(N_m-1):N_m;
    %Lesser Green function between atoms m(1) and m(2) obtained from
    %GL=i*Gr*f_L*Gamma_L*Ga+i*Gr*f_R*Gamma_R*Ga
    gl=1i*gr(j1+(m(1)-1)*2,j1)*F_left*gamma_left*gr(j1+2*(m(2)-1),j1)'+1i*gr(j1+(m(1)-1)*2,j2)*F_right*gamma_right*gr(j1+(m(2)-1)*2,j2)';
    %Current between atoms m(1) and m(2) obtained from real part of G_L 
    I(i)=-2*real(trace(S*gl));
end
end
function DS=dos_2d(E,t)
DS=zeros(size(E));
for i=1:max(size(E))
    if E(i)==0 || E(i)>4*t || E(i)<-4*t
        DS(i)=0;
    else
        DS(i)=1/(4*pi^2*t)*ellipke(sqrt(1-(E(i)/(4*t)).^2));%analytical formula for 2D DOS involves elliptic integral of the first kind
    end
end
end
function [out]=transmission(sigma_left,sigma_right,gr)
%computes transmission function for 1D tight-binding chain attached to ferromagnetic electrodes

N_m=length(gr);

gamma_left=1i*(sigma_left-sigma_left'); %level broadening matrix from left F electrode
gamma_right=1i*(sigma_right-sigma_right'); %level broadening matrix from right F electrode

ga=gr'; %advanced green function

out=real(trace(gamma_right*gr(N_m-1:N_m,1:2)*gamma_left*ga(1:2,N_m-1:N_m))); %Caroli et al. formula

end %function transmission 
function [Ich,ISx,ISy,ISz,T]=observables(d,V_bias,t,Tmpr,fermi_energy,hamil,sigma_left,sigma_right)

    mu=[fermi_energy,fermi_energy+V_bias];%chemical potential of the two leads. You can 
    %distribute V_bias between two leads anyway you want. Here we shifted the potential of right lead 

    if d==1
        T=Trans(fermi_energy,hamil,sigma_left,sigma_right);
        %Define a function that gets as input the energy of the electron and a 2x2 matrix for the
        %spin orientation of the calculated current: 
        %For charge current S=eye(2), 
        %for spin-x current S=[0,1;1,0], 
        %for spin-y current S=[0,-1i;1i,0], 
        %for spin-z current S=[1,0;0,-1], 
        F_I=@(Energy,S)arrayfun(@(E)Current(S,mu,Tmpr,E,0,hamil,sigma_left,sigma_right,[2,3]),Energy);
        %By uncommenting following part you can plot current vs energy of
        %the electrons. Notice that only electrons arround Fermi energy
        %contribute to the charge current.
%         x=(-6:0.01:6);
%         plot(x,F_I(x,eye(2)))
%         pause() %press any key to do the next calculation
        
    elseif d==2
        T=quadgk(@(x)Trans(fermi_energy-2*t*cos(x),hamil,sigma_left,sigma_right),-pi,pi)/(2*pi);
        F_I=@(Energy,S)arrayfun(@(E)quadgk(@(x)Current(S,mu,Tmpr,E,2*t*cos(x),hamil,sigma_left,sigma_right,[2,3])/(2*pi),-pi,pi),Energy);
%         x=(-6:0.01:6);
%         plot(x,F_I(x,eye(2)))
%         pause() %press any key to do the next calculation
    elseif d==3
        %uncomment following part if you want to use trapzoidal rule to do
        %the integration. Here we use adaptive Gauss-Kronrod quadrature
        %method to do the integration.
%         x=-4:de:4; %loop over in-plane kinetic energy grid
%         T_parallel=trapz(x,dos_2d(x,t).*Trans(fermi_energy-x,hamil,sigma_left,sigma_right));
        T=quadgk(@(x)dos_2d(x,t).*Trans(fermi_energy-x,hamil,sigma_left,sigma_right),-4.1,4);
        F_I=@(Energy,S)arrayfun(@(E)quadgk(@(x)dos_2d(x,t).*Current(S,mu,Tmpr,E,x,hamil,sigma_left,sigma_right,[2,3])/(2*pi),-4.1,4),Energy);
%         x=(-6:0.01:6);
%         plot(x,F_I(x,eye(2)))
%         pause() %press any key to do the next calculation
    end
    %Here we perform integration over energy to find total spin-z
    %current. Pay attention to the range of integration because in the
    %noncolinear case you have to integrate from the bottom of the band
    %(unlike here that only energies arrounf Fermi energy contribute to the result)
    if V_bias==0
        Ich=0;ISx=0;ISy=0;ISz=0;
    else
        Range=[fermi_energy-10*V_bias,fermi_energy+10*V_bias];
        Ich=quadgk(@(E)F_I(E,eye(2)),Range(1),Range(2));
        ISx=-quadgk(@(E)F_I(E,[0,1;1,0]),Range(1),Range(2));
        ISy=-quadgk(@(E)F_I(E,[0,-1i;1i,0]),Range(1),Range(2));
        ISz=-quadgk(@(E)F_I(E,[1,0;0,-1]),Range(1),Range(2));
    end
end

