%Program to set up Hamiltonian of finite-size graphene lattice
%with nearest neighbor hopping and possibly periodic boundary conditions

%%
clear all

%Parameters:
%%
%Device size
N_x=7; %Number of carbon atoms along the x-axis (should be odd for horizontal zigzag edge) 
N_y=4; %Number of carbon atoms along the y-axis (should be even for vertical armchair edge)
modes=N_y;

%%
%Hamiltonian parameters
gamma=1; %Unit of energy
t=-gamma; %hopping inside the sample
E0=0*gamma; %on-site potential


%%
Hsize=N_x*N_y; %Size of the full Hamiltonian matrix
H(Hsize,Hsize)=0; %Fill the full Hamiltonian matrix with zeros

%Hamiltonian of individual rows of atoms, such as:
% H_11    or    H_22
% C             C
% |
% C             C
%               |
% C             C
% |
% C             C

H_11=E0*diag(ones(1,N_y))+t*diag(ones(1,N_y-1),1)+t*diag(ones(1,N_y-1),-1); %First type of rows of atoms in the brick lattice
for n=2:2:N_y-1 %Remove hoppings in H_11 according to the brick lattice 
    H_11(n,n+1)=0;
    H_11(n+1,n)=0;
end
    
H_22=E0*diag(ones(1,N_y))+t*diag(ones(1,N_y-1),1)+t*diag(ones(1,N_y-1),-1); %Second type of rows in the brick lattice
for n=1:2:N_y %Remove hoppings in H_22 according to the brick lattice 
    H_22(n,n+1)=0;
    H_22(n+1,n)=0;
end

H_ij=t*diag(ones(1,N_y),0); %Matrices describing coupling between the two rows of atoms

%%
%Introduce periodic boundary conditions:
H_11(1,N_y)=t; H_11(N_y,1)=t;
H_22(1,N_y)=t; H_22(N_y,1)=t;

H(1:N_y,Hsize-N_y+1:Hsize)=H_ij;
H(Hsize-N_y+1:Hsize,1:N_y)=H_ij;
%%
for n=1:2:N_x
    H(1+(n-1)*modes:n*modes,1+(n-1)*modes:n*modes)=H_11; %Put Hamiltonians of individual layers as diagonal blocks of the full Hamiltonian
end   

for n=2:2:N_x
    H(1+(n-1)*modes:n*modes,1+(n-1)*modes:n*modes)=H_22; %Put Hamiltonians of individual layers as diagonal blocks of the full Hamiltonian
end   

for n=2:N_x
    H(1+(n-1)*modes:n*modes,1+(n-2)*modes:(n-1)*modes)=H_ij; %Put coupling matrices as off-diagonal blocks
    H(1+(n-2)*modes:(n-1)*modes,1+(n-1)*modes:n*modes)=H_ij'; %This is the block symmetric to the previous one which is 
    %its Hermitian conjugate to ensure the full Hamiltonian is hermitian matrix
end

%%
%Visualize Hamiltonian matrix to check if it is correct
%visual_graphene_H(N_x,N_y,H);


%%
%Find DOS using Green functions
tic, eta=i*10^(-1); %Small imaginary number needed to define retarded G^r and advanced G^a Green functions; 
j=1;
for E=-3.5:0.06:3.5 %energy loop
    gr=inv((E+eta)*eye(Hsize)-H);
    dos(j)=-1/pi*trace((gr-gr')/(2*i)); %Im part of a matrix G is Im G = (G-G')/2i
    j=j+1;
end %energy loop
toc

%%
%Plot DOS versus Energy using the same grid of energy values as in the for loop above
figure(2);
E=-3.5:0.06:3.5; %Energy values at which dos will be plotted
plot(E,dos);
xlabel('Energy/t'); 
ylabel('Total DOS');


