Python Problem with coding the band structure of graphene nanoribbons

[fatemehae6339]

Hi, I have a Matlab code for calculating the band structure of graphene nanoribbon which is working fine, but I wanted to convert it to python and I've done it. I guess I made a mistake somewhere because the plot is totally wrong. Could someone check it for me, please?

 

[Master1022]

Can you provide the equivalent (correct) Matlab code? It will help to have a reference for comparison.

 

[fatemehae6339]

Yeah, of course. Thank you in advance

 

[Master1022]

In the Python code, you have:

H_tot=H00+H01*exp(1j*kx*3*a)+H01.T*exp(-1j*kx*3*a)

whereas in Matlab you have:

 H=Ham+h01'*exp(1j*kx*3*a)+h01*exp(-1j*kx*3*a);

Does H01 = h01? If so, then I think you have used the wrong sign for the exponents ($\pm j$) in the exponentials?

 

[fatemehae6339]

You're right, I made a mistake there. I corrected that, but it didn't change the result

 

[Master1022]

Okay, the next step is to use the debugging tool within the Matlab and set break points at each point within that FOR loop. Does $d$ in Matlab come out exactly the same as $Ev$ in Python? If not, then compare (between Python and Matlab) to see whether each of the lines in the FOR loop output a different result. There are only 4 lines of code within that FOR loop so you should be able to determine which, if any, are causing the discrepancies.

If it isn't that, then it will likely be something to do with the plotting code as that is the only bit of code left in Python.

 

[fatemehae6339]

Obviously not, they're not the same. I think the problem is in the calculating of H_tot. the final matrix is different from the same matrix in Matlab. but I don't know how to correct it.
The procedure is,
First, I must calculate total Hamiltonian, i.e., H_tot for each kx. Then find the eigenvalues of each matrix and put them in each column of a matrix named Ev. Afterward, plot each row in terms of kx values.

and since I am new in python I don't know where I made a mistake

 

[Master1022]

I think the problem is in the calculating of H_tot. the final matrix is different from the same matrix in Matlab.

Which components are different? Did you inspect the individual components within that line to ensure that they are yielding the same values in both languages? There are two matrices (and a transpose) and two constants. Also, in your Python code you aren't iterating through the KX list. You can change the line to:

for kx in Kpoints:

I am only suggesting this as you have written a note in the Python that you are sure of all the code up until final portion. Therefore, if something has gone wrong it will be at the end.

[EDIT]: I used the Matlab name instead by mistake. Have changed KX to Kpoints
[EDIT #2]: Note, if you do this method, you do NOT need to change the H_tot line as kx will refer to the actual values within the elements. The method BELOW references the INDICES of elements in the array and thus requires H_tot to be changed

 

[Master1022]

Or if you want to use range in python, note that range(n) goes from 0 to $n - 1$ and indexing starts from 0 (unlike Matlab where it starts from 1). Therefore, you could also change the two lines to:

for i in range (len(Kpoints)):

(example, an array of 4 items will have a length 4 and the above will loop from 0 to 3 which will access all the indices of the elements in the array)

If you are going to use above, then you also need to change your H_tot line to:

 H_tot=H00+H01*exp(1j*Kpoints[i]*3*a)+H01.T*exp(-1j*Kpoints[i]*3*a)

 

[fatemehae6339]

Which components are different? Did you inspect the individual components within that line to ensure that they are yielding the same values in both languages? There are two matrices (and a transpose) and two constants. Also, in your Python code you aren't iterating through the KX list. You can change the line to:

for kx in Kpoints:

I am only suggesting this as you have written a note in the Python that you are sure of all the code up until final portion. Therefore, if something has gone wrong it will be at the end.

[EDIT]: I used the Matlab name instead by mistake. Have changed KX to Kpoints
[EDIT #2]: Note, if you do this method, you do NOT need to change the H_tot line as kx will refer to the actual values within the elements. The method BELOW references the INDICES of elements in the array and thus requires H_tot to be changed

 

[fatemehae6339]

Hello again, and thank you very much
that day I didn't really understand what you mean by kx=Kx. I understood yesterday and changed my code this way:

for kx in Kpoints:
H_tot=H00+H01*exp(-1j*kx*3*a)+H01.T*exp(1j*kx*3*a)
ii=0
Ev[:,ii]=np.linalg.eigvals(H_tot)
Ev=np.sort(Ev)
ii=ii+1

but the result still incorrect

 

[fatemehae6339]

Or if you want to use range in python, note that range(n) goes from 0 to $n - 1$ and indexing starts from 0 (unlike Matlab where it starts from 1). Therefore, you could also change the two lines to:

for i in range (len(Kpoints)):

(example, an array of 4 items will have a length 4 and the above will loop from 0 to 3 which will access all the indices of the elements in the array)

If you are going to use above, then you also need to change your H_tot line to:

 H_tot=H00+H01*exp(1j*Kpoints[i]*3*a)+H01.T*exp(-1j*Kpoints[i]*3*a)

I tried this way,too:

for i in range(len(Kpoints)):

H_tot=H00+H01*exp(-1j*Kpoints*3*a)+H01.T*exp(1j*Kpoints*3*a)
Ev[:,i]=np.linalg.eigvals(H_tot)
Ev=np.sort(Ev)

and still the same result

 

[Master1022]

Hi,

I have responded to the code in both methods - feel free to choose either one.

Reference to method #1

for kx in Kpoints:
H_tot=H00+H01*exp(-1j*kx*3*a)+H01.T*exp(1j*kx*3*a)
ii=0
Ev[:,ii]=np.linalg.eigvals(H_tot)
Ev=np.sort(Ev)
ii=ii+1

but the result still incorrect

You are accessing the correct values of kx now. This method basically iterates through the elements (i.e. not the indices) of the list. However, within the loop, you have set

 ii = 0

during each iteration and therefore you will always be accessing Ev[:,0]. I think that definition ought to go outside the loop definition:

ii=0
for kx in Kpoints:
    H_tot=H00+H01*exp(-1j*kx*3*a)+H01.T*exp(1j*kx*3*a)
    Ev[:,ii]=np.linalg.eigvals(H_tot)
    Ev=np.sort(Ev)
    ii=ii+1

Reference to Method #2

I tried this way,too:

for i in range(len(Kpoints)):

H_tot=H00+H01*exp(-1j*Kpoints*3*a)+H01.T*exp(1j*Kpoints*3*a)
Ev[:,i]=np.linalg.eigvals(H_tot)
Ev=np.sort(Ev)

and still the same result

This should instead read:

for i in range(len(Kpoints)):
    H_tot=H00+H01*exp(-1j*Kpoints[i]*3*a)+H01.T*exp(1j*Kpoints[i]*3*a)
    Ev[:,i]=np.linalg.eigvals(H_tot)
    Ev=np.sort(Ev)

So in this method, we are basically iterating through the indices of the list and store the index in the variable i. Therefore, to access the $i^{th}$ element, we need to do

Kpoints[i]

If that still doesn't fix it, it is a good idea to look at each element of the H_tot line and look at what values each variable is taking. One of them will be different that the corresponding values in Matlab and it should help to point the errors out from there. Unfortunately, I don't know much about the physics here so I can't check the equations, but I do know some Python/Matlab. Also, do you need to np.sort the matrix Ev within each iteration of the for loop (as opposed to after the loop)?

Let me know if that didn't make sense.

 

[fatemehae6339]

Actually, when I commet this line(np.sort...), the result is much more closer to the correct answer, but still incorrect.

 

[Master1022]

Actually, when I commet this line(np.sort...), the result is much more closer to the correct answer, but still incorrect.

Hi,

Is the H_tot line now outputting the same value in both Python and Matlab? If so, then what about the line that grabs the eigenvalues? I think isolating where the error lies in a sequential manner is probably the best thing to do.

In terms of using the sort function, just to confirm: do you want to sort by the real parts of the eigenvalues in ascending order?

 

[fatemehae6339]

Hi,

Thank you very very much for your help.
I add axis option to the sort command:

Ev=np.sort(Ev,axis=0)

and, the result is now correct.