Graphene Sparse vs Dense  Energy calculation values dependent on graphene length
This is related to my other post, but as I see them as different (maybe not?) problems, I will post a second thread.
The situation is, again, a graphene strip with periodic boundary conditions in one direction and open boundaries in the other.
Using the sparse solver from linalg, I have solved for energies and gotten reasonable results for the energies, except for the case when the length (along the axis of the open boundary condition) is metallic. Attached is a code that shows it for length = 17, which is a metallic length for graphene (has energies at zero), but it does not produce messy energy calculation for when it is not a metallic length, i.e. something like length = 16.
""" Graphene wire with periodic boundary conditions in the vertical direction. """
import kwant from math import pi, sqrt, tanh, cos, ceil, floor, atan, acos, asin from cmath import exp import numpy as np import scipy import scipy.linalg as lina
sin_30, cos_30, tan_30 = (1 / 2, sqrt(3) / 2, 1 / sqrt(3))
def create_closed_system(length, width, lattice_spacing, onsite_potential, hopping_parameter, boundary_hopping):
padding = 0.5*lattice_spacing*tan_30
def calc_total_length(length): total_length = length N = total_length//lattice_spacing # Number of times a graphene hexagon fits in horizontially fully new_length = N*lattice_spacing + lattice_spacing*0.5 diff = total_length  new_length if diff != 0: length = length  diff total_length = new_length return total_length
def calc_width(width,lattice_spacing, padding): stacking_width = lattice_spacing*((tan_30/2)+(1/(2*cos_30))) N = width//stacking_width if N % 2 == 0.0: # Making sure that N is odd. N = N1 new_width = N*stacking_width + padding width = new_width return width, int(N)
def rectangle(pos): x,y = pos if (0 < x <= total_length) and (padding <= y <= width padding): return True return False
def lead_shape(pos): x, y = pos return 0  padding <= y <= width
def tag_site_calc(x): return int(1*(x*0.5+0.5))
#Initation of geometrical limits of the lattice of the system total_length = calc_total_length(length) width, N = calc_width(width, lattice_spacing, padding)
# The definition of the potential over the entire system def potential_e(site,pot): return pot
### Definig the lattices ###
graphene_e = kwant.lattice.honeycomb(a=lattice_spacing,name='e') a, b = graphene_e.sublattices
sys = kwant.Builder()
# The following functions are required for input in the def onsite_shift_e(site, pot): return potential_e(site,pot)
sys[graphene_e.shape(rectangle, (0.5*lattice_spacing, 0))] = onsite_shift_e sys[graphene_e.neighbors()] = hopping_parameter
### Boundary conditions for scattering region ### for site in sys.sites(): (x,y) = site.tag if float(site.pos[1]) < 0: if str(site.family) == "<Monatomic lattice e1>": sys[b(x,y),a(int(x+tag_site_calc(N)),N)] = boundary_hopping
kwant.plot(sys)
return sys
def eigen_vectors_and_values(sys,sparse_dense,k): #sparse = 0, dense = 1 if sparse_dense == 0: ham_mat = sys.hamiltonian_submatrix(params=dict(pot=0.0), sparse=True) eigen_val, eigen_vec = scipy.sparse.linalg.eigsh(ham_mat.tocsc(), k=k, sigma=0, return_eigenvectors=True) if sparse_dense == 1: ham_mat = sys.hamiltonian_submatrix(params=dict(pot=0.0), sparse=False) eigen_val, eigen_vec = lina.eigh(ham_mat)
#sort the ee and ev in ascending order idx = eigen_val.argsort() eigen_val= eigen_val[idx] eigen_vec = eigen_vec[:,idx]
if sparse_dense == 1: eigen_val = eigen_val[int(len(eigen_val)*0.5k*0.5):] # Remove first part eigen_val = eigen_val[:k] # Take first kvalues eigen_vec = eigen_vec[:,int(len(eigen_val)*0.5k*0.5):] eigen_vec = eigen_vec[:,:k]
return eigen_vec, eigen_val
def main():
sys = create_closed_system(length = 17.0, width = 20.0, lattice_spacing = 1.0, onsite_potential = 0.0, hopping_parameter = 1.0, boundary_hopping = 1.0)
sys = sys.finalized()
#sparse = 0, dense = 1 eigen_vec, eigen_val = eigen_vectors_and_values(sys,sparse_dense=0,k=50) eigen_vec2, eigen_val2 = eigen_vectors_and_values(sys,sparse_dense=1,k=50)
print("Eigenenergy difference (between sparse and dense): ") for i_ in range(len(eigen_val)): print(eigen_val[i_]  eigen_val2[i_] )
return
if __name__ == "__main__": main()
Dear Henry,
Your code is full of parts with potential mistakes. I mean by that, using scipy.sparse.linalg.*eigsh *to get the eigenvalues and eigenvectors might be a source of error because you only calculate part of the eigenvalues, and then need to reorder them in the same way linalg.eigh computes all the eigenvalues. What you should notice is that you have many degenerate eigenvalues and reordering them will not give you the requested order for the eigenvectors. Another point to add concerning the degenerate eigenvalues: the eigenvector is any combination of linearly independent eigenvectors which means it is not unique.
The comparison between the eigenvectors is not safe. (even a simple global phase can make it not safe!)
I hope this helps , Adel
On Tue, Sep 28, 2021 at 3:07 PM Henry Axt henry.axt@gmail.com wrote:
This is related to my other post, but as I see them as different (maybe not?) problems, I will post a second thread.
The situation is, again, a graphene strip with periodic boundary conditions in one direction and open boundaries in the other.
Using the sparse solver from linalg, I have solved for energies and gotten reasonable results for the energies, except for the case when the length (along the axis of the open boundary condition) is metallic. Attached is a code that shows it for length = 17, which is a metallic length for graphene (has energies at zero), but it does not produce messy energy calculation for when it is not a metallic length, i.e. something like length = 16.
""" Graphene wire with periodic boundary conditions in the vertical direction. """
import kwant from math import pi, sqrt, tanh, cos, ceil, floor, atan, acos, asin from cmath import exp import numpy as np import scipy import scipy.linalg as lina
sin_30, cos_30, tan_30 = (1 / 2, sqrt(3) / 2, 1 / sqrt(3))
def create_closed_system(length, width, lattice_spacing, onsite_potential, hopping_parameter, boundary_hopping):
padding = 0.5*lattice_spacing*tan_30 def calc_total_length(length): total_length = length N = total_length//lattice_spacing # Number of times a graphene
hexagon fits in horizontially fully new_length = N*lattice_spacing + lattice_spacing*0.5 diff = total_length  new_length if diff != 0: length = length  diff total_length = new_length return total_length
def calc_width(width,lattice_spacing, padding): stacking_width = lattice_spacing*((tan_30/2)+(1/(2*cos_30))) N = width//stacking_width if N % 2 == 0.0: # Making sure that N is odd. N = N1 new_width = N*stacking_width + padding width = new_width return width, int(N) def rectangle(pos): x,y = pos if (0 < x <= total_length) and (padding <= y <= width padding): return True return False def lead_shape(pos): x, y = pos return 0  padding <= y <= width def tag_site_calc(x): return int(1*(x*0.5+0.5)) #Initation of geometrical limits of the lattice of the system total_length = calc_total_length(length) width, N = calc_width(width, lattice_spacing, padding) # The definition of the potential over the entire system def potential_e(site,pot): return pot ### Definig the lattices ### graphene_e = kwant.lattice.honeycomb(a=lattice_spacing,name='e') a, b = graphene_e.sublattices sys = kwant.Builder() # The following functions are required for input in the def onsite_shift_e(site, pot): return potential_e(site,pot) sys[graphene_e.shape(rectangle, (0.5*lattice_spacing, 0))] =
onsite_shift_e sys[graphene_e.neighbors()] = hopping_parameter
### Boundary conditions for scattering region ### for site in sys.sites(): (x,y) = site.tag if float(site.pos[1]) < 0: if str(site.family) == "<Monatomic lattice e1>": sys[b(x,y),a(int(x+tag_site_calc(N)),N)] =
boundary_hopping
kwant.plot(sys) return sys
def eigen_vectors_and_values(sys,sparse_dense,k): #sparse = 0, dense = 1 if sparse_dense == 0: ham_mat = sys.hamiltonian_submatrix(params=dict(pot=0.0), sparse=True) eigen_val, eigen_vec = scipy.sparse.linalg.eigsh(ham_mat.tocsc(), k=k, sigma=0, return_eigenvectors=True) if sparse_dense == 1: ham_mat = sys.hamiltonian_submatrix(params=dict(pot=0.0), sparse=False) eigen_val, eigen_vec = lina.eigh(ham_mat)
#sort the ee and ev in ascending order idx = eigen_val.argsort() eigen_val= eigen_val[idx] eigen_vec = eigen_vec[:,idx] if sparse_dense == 1: eigen_val = eigen_val[int(len(eigen_val)*0.5k*0.5):] # Remove
first part eigen_val = eigen_val[:k] # Take first kvalues eigen_vec = eigen_vec[:,int(len(eigen_val)*0.5k*0.5):] eigen_vec = eigen_vec[:,:k]
return eigen_vec, eigen_val
def main():
sys = create_closed_system(length = 17.0, width = 20.0, lattice_spacing = 1.0, onsite_potential = 0.0, hopping_parameter = 1.0, boundary_hopping = 1.0) sys = sys.finalized() #sparse = 0, dense = 1 eigen_vec, eigen_val =
eigen_vectors_and_values(sys,sparse_dense=0,k=50) eigen_vec2, eigen_val2 = eigen_vectors_and_values(sys,sparse_dense=1,k=50)
print("Eigenenergy difference (between sparse and dense): ") for i_ in range(len(eigen_val)): print(eigen_val[i_]  eigen_val2[i_] ) return
if __name__ == "__main__": main()
Dear Abbout,
Thank you for the response!
I understand what you mean by the reordering of the eigenvalues/eigenstates, but I believe I deal with this problem with the following lines:
idx = eigen_val.argsort() eigen_val= eigen_val[idx] eigen_vec = eigen_vec[:,idx]
This makes me reaslise that the comparison of the eigenvectors calculated by sparse and dense is not sensible. But it seems that the calculation of the eigenvalues is sensible and the length dependent success of the sparse algorithm in the calculation of the eigenvalues still seems very mysterious to me.
participants (2)

Abbout Adel

Henry Axt