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()

Reply via email to