% spinNEGF1D
% Rashba term: sx*ky-sy*kx = (1/2/i/a) (e*ikya-e*-ikya)*sx-(e*ikxa-e*-ikxa)*sy
clear all
sx=[0 1;1 0];sy=[0 -i;i 0];sz=[1 0;0 -1];
[vx,dd]=eig(sx);[vy,dd]=eig(sy);[vz,dd]=eig(sz);
V=vx;dn=V*[1 0;0 0]*V';up=V*[0 0;0 1]*V';
 
% Inputs
hbar=1.06e-34;q=1.6e-19;m=0.03*9.1e-31;
a=1e-9;t0=(hbar^2)/(2*m*(a^2)*q);zplus=i*1e-12;
Np=151;L=zeros(Np);R=L;L(1,1)=1;R(Np,Np)=1;X=a*[0:1:Np-1];
    M=1*kron(eye(Np),[1 1]);C=eye(2*Np);
    Sx=kron(eye(Np),sx);Sy=kron(eye(Np),sy);Sz=kron(eye(Np),sz);
 
EE=0.2;rashba=0e-11;Bfield=3e-2;
    alpha=2*t0*eye(2);beta=-t0*eye(2)-(rashba/i/2/a)*sy;
    N3=ceil(Np/3);
    bar=[zeros(1,N3),ones(1,N3),zeros(1,Np-2*N3)];
    H0=kron(eye(Np),alpha)+kron(diag(bar),Bfield*sz);
        Hu=kron(diag(ones(1,Np-1),+1),beta);
            Hd=kron(diag(ones(1,Np-1),-1),beta');
    H=H0+Hd+Hu;
    
    ck=(1-(EE+zplus)/(2*t0));ka=acos(ck);
    sigu=(-t0*exp(i*ka))*up;sigd=(-t0*exp(i*ka))*dn;
    sig1u=kron(L,sigu);sig1d=kron(L,sigd);  
    sig2u=kron(R,sigu);sig2d=kron(R,sigd);  
    gam1u=i*(sig1u-sig1u');gam1d=i*(sig1d-sig1d');
    gam2u=i*(sig2u-sig2u');gam2d=i*(sig2d-sig2d');
    
    G=inv((EE*eye(2*Np))-H-sig1u-sig1d-sig2u-sig2d);
    A=i*(G-G');Gn=G*gam1u*G';
    
    DGn=Gn.*kron(eye(Np),ones(2));
    n=M*real(diag(DGn*C));nx=M*real(diag(DGn*Sx));
    ny=M*real(diag(DGn*Sy));nz=M*real(diag(DGn*Sz));
 
    DJ=-i*(Hd*Gn-Gn*Hu).*kron(eye(Np),ones(2));
    DJR=i*(Hu*Gn-Gn*Hd).*kron(eye(Np),ones(2));
    DJC=i*(H0*Gn-Gn*H0).*kron(eye(Np),ones(2));
    J=M*real(diag(DJ*C));
    Jx=M*real(diag(DJ*Sx));JxR=M*real(diag(DJR*Sx));JxC=M*real(diag(DJC*Sx));
    Jy=M*real(diag(DJ*Sy));JyR=M*real(diag(DJR*Sy));JyC=M*real(diag(DJC*Sy));
    Jz=M*real(diag(DJ*Sz));JzR=M*real(diag(DJR*Sz));JzC=M*real(diag(DJC*Sz));
 
    % Checking: order should not matter
    n1=M*real(diag(C*DGn));nx1=M*real(diag(Sx*DGn));
    ny1=M*real(diag(Sy*DGn));nz1=M*real(diag(Sz*DGn));
    J1=M*real(diag(C*DJ));Jx1=M*real(diag(Sx*DJ));
    Jy1=M*real(diag(Sy*DJ));Jz1=M*real(diag(Sz*DJ));

hold on
%% Spin current Jx (left/ on-site/ right)
figure(1)
h=plot(X,JxC,'m');hold on
set(h,'linewidth',[3.0])
h=plot(X,Jx,'rx');hold on
h=plot(X,JxR,'b');
set(h,'linewidth',[3.0])
set(gca,'Fontsize',[36])
grid on
% ylim([-0.1 0.1])
xlabel('distance[m]');
ylabel('Spin current J_x');
torque_x=sum(JxC)
%% Spin current Jy (left/ on-site/ right)
figure(2)
h=plot(X,JyC,'m');hold on
set(h,'linewidth',[3.0])
h=plot(X,Jy,'rx');hold on
h=plot(X,JyR,'b');
set(h,'linewidth',[3.0])
set(gca,'Fontsize',[36])
grid on
xlabel('distance[m]');
ylabel('Spin current J_y');
torque_y=sum(JyC)
%% Spin current Jz (left/ on-site/ right)
figure(3)
h=plot(X,JzC,'m');hold on
set(h,'linewidth',[3.0])
h=plot(X,Jz,'rx');hold on
h=plot(X,JzR,'b');
set(h,'linewidth',[3.0])
set(gca,'Fontsize',[36])
grid on
xlabel('distance[m]');
ylabel('Spin current J_z');
torque_z=sum(JzC)
%% Spin current Jx/Jy/Jz
P=J;P1=J1;Px=Jx;Px1=Jx1;Py=Jy;Py1=Jy1;Pz=Jz;Pz1=Jz1;
figure(4)
% h=plot(X,P,'y');hold on
% set(h,'linewidth',[3.0])
h=plot(X,P1,'co');hold on
% h=plot(X,Px,'r');
set(h,'linewidth',[3.0])
h=plot(X,Px1,'rx');
% h=plot(X,Py,'b');
set(h,'linewidth',[3.0])
h=plot(X,Py1,'b+');
% h=plot(X,Pz,'m');
set(h,'linewidth',[3.0])
h=plot(X,Pz1,'mo');
set(gca,'Fontsize',[36])
xlabel('distance[m]');
ylabel('Spin current J_x,_y,_z');
grid on
%% Spin density Sx/Sy/Sz
P=n;P1=n1;Px=nx;Px1=nx1;Py=ny;Py1=ny1;Pz=nz;Pz1=nz1;
figure(5)
% h=plot(X,P,'y');hold on
% set(h,'linewidth',[3.0])
h=plot(X,P1,'co');hold on
% h=plot(X,Px,'r');
set(h,'linewidth',[3.0])
h=plot(X,Px1,'rx');
% h=plot(X,Py,'b');
set(h,'linewidth',[3.0])
h=plot(X,Py1,'b+');
% h=plot(X,Pz,'m');
set(h,'linewidth',[3.0])
h=plot(X,Pz1,'mo');
set(gca,'Fontsize',[36])
xlabel('distance[m]');
ylabel('Spin density S_x,_y,_z');
grid on