Dear Adel,
Thank you very much for pointing out the thread
https://www.mail-archive.com/kwant-discuss@python.org/msg00278.html, I
completely missed it.
I studied the Silicene code by Koustav that had the same problem of my point
1): RuntimeError("Numbers of left- and right-propagating " RuntimeError:
Numbers of left- and right-propagating modes differ, possibly due to a
numerical instability.
Despite Koustav replied that the proposed solution works, I find that his code
(plus the sym.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
corrections) does not give consistent results. In particular KWANT cannot
calculate the local density of states for every choice of the dimensions of the
system and for every energy. The error it gives is the same, i.e. the numbers
of left- and right-propagating modes differ.
I also tried implementig the add_site_family in my calculation and, as in the
Silicene code by Koustav, running the code gives inconsistent results in
evaluating the wave function.
My conclusion is that the suggested solution might work for certain cases but
it does not scale well with the dimensions of the system.
There should be another way to overcome this problem, but I couldn't find it.
For reference and convenience, I attach the code I'm using at the end of this
message.
Thank you again for you help!
All the best,
Lorenzo
##############################################################
import kwant
import numpy as np
import math
from tqdm import tqdm
import matplotlib.pyplot as plt
import plotly
from mpl_toolkits.mplot3d import Axes3D
from boiadeh import *
from matplotlib import colormaps
from matplotlib.colors import Normalize
#Definition of the primitive vectors
lattice_constant= 1#0.246 #2.46 \Angstrom
armchair = False
zigzag = False
FF = True
if armchair:
a1, a2, d1 = [
lattice_constant*np.array([np.sqrt(3)/2,1./2]),
lattice_constant*np.array([np.sqrt(3)/2,-1./2]),
(lattice_constant/np.sqrt(3))*np.array([1./2,-np.sqrt(3)/2])]
elif zigzag:
#ZIGZAG
a1, a2, d1 = [
lattice_constant*np.array([1.,0]),
lattice_constant*np.array([1/2, np.sqrt(3)/2]),
(lattice_constant/np.sqrt(3))*np.array([0,1])]
elif FF:
a1, a2, d1 = [
lattice_constant*np.array([np.sqrt(3)/2,1./2]),
lattice_constant*np.array([np.sqrt(3)/2,-1./2]),
(lattice_constant/np.sqrt(3))*np.array([1,0])]
# Define the graphene lattice
graphene = kwant.lattice.general([a1,a2],[(0,0), d1], norbs=1)
A1, B1 = graphene.sublattices
sublattices = {
(1,0): A1,
(1,1): B1
}
nnn_hoppings_a = (((-1, 0), A1, A1), ((0, 1), A1, A1), ((1, -1), A1, A1))
nnn_hoppings_b = (((1, 0), B1, B1), ((0, -1), B1, B1), ((-1, 1), B1, B1))
sys = kwant.Builder()
#On site and hopping energies
on_site_sc = 0
on_site_lead = 0.
t=10#2.7
print('On site and hopping energies: \n')
print('Scatterer on site energy: ' + str(on_site_sc) + ' \n')
print('Leads on site energy: ' + str(on_site_lead) + ' \n')
print('NN hopping graphene: ' + str(t) + ' \n')
#Intra layer parameters SOC
lambda0= 0.1 # SOC A1-B1 intralayer
#Print selected SOC parameters
print('SOC Hamiltonian parameters:\n')
print('lambda_0 = ' + str(lambda0) +'\n')
#Print selected legth parameters
print('Length scales parameters:\n')
print('lattice constant = ' + str(lattice_constant) +'\n')
L=30*lattice_constant
W=20*lattice_constant#/np.sqrt(3)
#Definition of the shape
def rectangle(pos):
x, y = pos
return np.abs(x) <= L/2 and np.abs(y) <= W/2
#Define the sites of the ribbon
on_site_sc = 0.1*np.sqrt(3)*3*(0)*t
def onsite(site):
if site.family == A1:
return on_site_sc
elif site.family == B1:
return - on_site_sc
sys[graphene.shape(rectangle, (0, 0))]= onsite #on-site energy, fermi level
def family_colors(site):
if site.family == A1:
return 'r'
elif site.family == B1:
return 'g'
#nearest-neighbor hopping
sys[graphene.neighbors(1)] = t/2 #+ 0.5*1j*lambda0*(sigma_x-1j*sigma_y)
sys.eradicate_dangling()
t2=0.1*t
phi = -1*np.pi/4
#second nearest-neighbor hopping
sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] =
t2/2*np.cos(phi) + 1j*t2/2*np.sin(phi)#np.exp(-1j*phi)
sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] =
t2/2*np.cos(phi) + 1j*t2/2*np.sin(phi)#np.exp(1j*phi)
#Plot the system with attached leads
kwant.plot(sys,
site_size=0.18,
site_lw=0.01,
hop_lw=0.05,
site_color=family_colors,
fig_size = (8,8),
dpi=300
);
def lead_shape(pos):
x, y = pos
return np.abs(x) <= L/2
#### Define the leads ####
sym_left = kwant.TranslationalSymmetry(graphene.vec((1,-1)))
sym_right = kwant.TranslationalSymmetry(graphene.vec((-1,1)))
sym_left.add_site_family(graphene.sublattices[0],other_vectors=[(1, 1)])
sym_left.add_site_family(graphene.sublattices[1],other_vectors=[(1, 1)])
sym_right.add_site_family(graphene.sublattices[0],other_vectors=[(1, 1)])
sym_right.add_site_family(graphene.sublattices[1],other_vectors=[(1, 1)])
left_lead = kwant.Builder(sym_left)#, conservation_law=-sigma_z)
right_lead = kwant.Builder(sym_right)#, conservation_law=-sigma_z)
left_lead[graphene.shape(lead_shape, (0, 0))] = onsite # 4 * t*sigma_0
right_lead[graphene.shape(lead_shape, (0, 0))] = onsite # 4 * t*sigma_0
left_lead[graphene.neighbors(1)] = +t/2 #+ 0.5*1j*lambda0*(sigma_x-1j*sigma_y)
right_lead[graphene.neighbors(1)] = +t/2 #+ 0.5*1j*lambda0*(sigma_x-1j*sigma_y)
left_lead.eradicate_dangling()
right_lead.eradicate_dangling()
left_lead[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]]
= +t2/2*np.cos(phi) + 1j*t2/2*np.sin(phi)#1j*phi#np.exp(-1j*phi)
left_lead[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]]
= +t2/2*np.cos(phi) + 1j*t2/2*np.sin(phi)#1j*phi#np.exp(1j*phi)
right_lead[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]]
= +t2/2*np.cos(phi) + 1j*t2/2*np.sin(phi)#1j*phi#np.exp(-1j*phi)
right_lead[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]]
= +t2/2*np.cos(phi) + 1j*t2/2*np.sin(phi)#1j*phi#np.exp(1j*phi)
#Attach every lead
sys.attach_lead(left_lead);
sys.attach_lead(right_lead);
#System finalization
syst = sys.finalized()
#Plot the system with attached leads
kwant.plot(sys,
site_size=0.18,
site_lw=0.01,
hop_lw=0.05,
site_color=family_colors,
fig_size = (8,8),
dpi=300
);
#Extract the wave function at fixed energy
wf_en = kwant.wave_function(syst, energy=-3)
wf_en(0).shape
##############################################################
