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