Trying to replicate Local density of states

I am trying to replicate Fig. 3 (b). I am getting figure 3(a) but not getting fig 3(b). I have attached my code here. Please someone help me to solve the issue.

For plotting I am following equation 7 and 8. The equations are attached here.

The whole paper is attached here for any clarification.

clc;
clear all;
close all;
t = 1.03;% in eV (Electron Volts)
a = 3.82;%in Angstroms
ka = linspace(-pi,pi,1000);
k = linspace(-pi,pi,1000)/a;
% Zig Zag Graphene N = 4
A = zeros(8,8);
A(1,2) = t;
A(2,1) = t;A(2,3) = t;
A(3,2) = t;A(3,4) = t;
A(4,3) = t;A(4,5) = t;
A(5,4) = t;A(5,6) = t;
A(6,5) = t;A(6,7) = t;
A(7,6) = t;A(7,8) = t;
A(8,7) = t;

% Defining the Connectivity between the Unit Cells Matrix
B = zeros(8,8);
B(2,1) =t;
B(3,4) = t;
B(6,5) =t;
B(7,8) =t;
for k = 1:length(ka)
[V1,D1] = (eig(A + B*exp(1i.*ka(k))+ B'.*exp(-1i.*ka(k))));
eigE1(k,:) = diag(D1);
hold on
end

figure(1)
plot(ka,eigE1,ka,-eigE1,'linewidth',1.5)
xlabel('ka');
ylabel('Energy in eV');
title('Energy Vs Ka for N = 4 Z-SNR');
box on
grid on
set(gca,'XTick',[-ka(length(ka)),0,ka(length(ka))]);
set(gca,'XTickLabel',{'-\pi','0','\pi'});
line([ka(length(ka)) ka(length(ka))],[4 -4]);
line([-ka(length(ka)) -ka(length(ka))],[4 -4]);

% Density of states
% Define the size of the matrix
rows = 1000;
cols = 8;
% Create a matrix with zeros
E = zeros(rows, cols);
% Define the diagonal elements
diagonal_values = linspace(-3, 3, min(rows, cols));
% Assign the diagonal values to the matrix
for i = 1:length(diagonal_values)
    E(i, i) = diagonal_values(i);
end
    gam=0.05;
    m=1;
    syms m
    DE1=symsum((2*gam)./((E-eigE1).^2+(gam).^2),m);
    i=i+1;
    %syms j
   % sum_DOS=symsum(1/(2*pi)*DE1,j,1,1);
   sum_DOS=(1/(2*pi)*sum(DE1,k));
    

hold off

figure(2)
plot(sum_DOS,E,'g');
xlabel('D(E)------>')
ylabel('E(eV)------>')
set(gca,'XLIM',[0,4])
box on;grid on

Answer 1

[Avni Agrawal]

Hi Swastik,

I understand that you want to density of states. The error message for second plot indicates that plot is being passed data that it cannot interpret as valid plotting data. Given that ‘sum_DOS’ is described as 1000*8 sym, it seems that ‘sum_DOS’is a symbolic matrix, not a numeric array. The plot function cannot handle symbolic expressions directly.

To resolve this issue, you'll need to convert the symbolic expressions to numeric values before plotting. You'll need to use the ‘subs’ function to replace all the symbolic variables in‘sum_DOS’ with their corresponding numerical values.

Before line 68 in your function, you can insert the following line to verify and modify the data accordingly:

figure(2)
% update sum_DOS using subs function for numerical values
sum_DOS = subs(sum_DOS, m, 1); 
plot(sum_DOS,E,'g');

Please refer to the following documentation page for more information on the ‘subs’ function: https://www.mathworks.com/help/symbolic/subs.html

I hope this helps.

[SWASTIK SAHOO]

Thank you so much. It really helps.