Hello, I am not too sure where my problem lies so, I will post my question here, since I might suspect it has something to do with the kwant interface. I have a simple graphene strip with periodic boundaries in the vertical direction and open boundaries in the other direction. To calculate the energies and states of this system, I feed the hamiltonian from kwant into either the dense or the sparse solver from linalg. I do not understand why I am getting significant difference in the values of the eigenstates using sparse vs using density, especially considering the small size of my system. Attached is a code that reproduces the problem: """ 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 = N-1 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.5-k*0.5):] # Remove first part eigen_val = eigen_val[:k] # Take first k-values eigen_vec = eigen_vec[:,int(len(eigen_val)*0.5-k*0.5):] eigen_vec = eigen_vec[:,:k] return eigen_vec, eigen_val def main(): sys = create_closed_system(length = 16.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: ") for i_ in range(len(eigen_val)): print(eigen_val[i_] - eigen_val2[i_] ) new_m = eigen_vec - eigen_vec2 new_m = abs(new_m) print("Max difference in eigenvector values: ", np.amax(new_m)) return new_m if __name__ == "__main__": new_m = main()