Dear Christoph, I simplified the posted program in order to track the error. It seems that it comes from this part: ###################################### for site in lead2.sites(): (x,y) = site.tag if str(site.family) == "<Monatomic lattice e0>": lead2[c(x,y),a(x,y)] =1 if str(site.family) == "<Monatomic lattice e1>": lead2[d(x,y),b(x,y)] = 1 ###################################### I could not go deeper to find the source of the error. The simplified program which shows the difference between"norbs=None" and "norbs=1" is pasted below: import kwant from matplotlib import pyplot from math import pi, sqrt, tanh, cos, ceil from cmath import exp import numpy as np import time import textwrap as tw sin_30, cos_30, tan_30 = (1 / 2, sqrt(3) / 2, 1 / sqrt(3)) def create_system(N=31, Norbs=None ): def rectangle(pos): x,y = pos if (0 < x <= 0.4*6.5) and (-3<= y <= 3): return True return False def lead_shape(pos): x, y = pos return - 3 <= y <= 3 graphene_e = kwant.lattice.honeycomb(a=1,name='e', norbs=Norbs) graphene_h = kwant.lattice.honeycomb(a=1,name='h', norbs=Norbs) a, b = graphene_e.sublattices c, d = graphene_h.sublattices sys = kwant.Builder() sys[graphene_e.shape(rectangle, (0.5, 0))] = 0.3 sys[graphene_h.shape(rectangle, (0.5, 0))] = -0.3 sys[graphene_e.neighbors()] = -1 sys[graphene_h.neighbors()] = 1 # Symmetry of the left and right leads sym_l = kwant.TranslationalSymmetry((-1, 0)) sym_r = kwant.TranslationalSymmetry((1, 0)) lead0 = kwant.Builder(sym_l) # Left electron lead0[graphene_e.shape(lead_shape, (1,1))] = -1 lead0[graphene_e.neighbors()] = -1 lead1 = kwant.Builder(sym_l) # Left hole lead1[graphene_h.shape(lead_shape, (1,1))] = 1 lead1[graphene_h.neighbors()] = 1 lead3 = kwant.Builder(sym_r) # Right electron lead3[graphene_e.shape(lead_shape, (1,1))] = -0.3 lead3[graphene_e.neighbors()] = -1 lead4 = kwant.Builder(sym_r) # Right hole lead4[graphene_h.shape(lead_shape, (1,1))] = 0.3 lead4[graphene_h.neighbors()] = 1 lead2 = kwant.Builder(sym_r) # Right superconducting lead2.update(lead3) lead2.update(lead4) # Hopping between electron and hole lead lattices in superconducting lead for site in lead2.sites(): (x,y) = site.tag if str(site.family) == "<Monatomic lattice e0>": lead2[c(x,y),a(x,y)] =1 if str(site.family) == "<Monatomic lattice e1>": lead2[d(x,y),b(x,y)] = 1 # Attaching the left electron, left hole and right superconducting lead sys.attach_lead(lead0) sys.attach_lead(lead1) sys.attach_lead(lead2) return sys, [lead0, lead1, lead2] def plot_conductance(syst, energies, lead_in, lead_out): data = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) data.append(smatrix.transmission(lead_out, lead_in)) pyplot.figure() pyplot.plot(energies, data) return data def main(): #case 1 sys, leads = create_system(Norbs=None ) sys = sys.finalized() kwant.plot(sys,fig_size=(20,20)) Es = np.linspace(0.0001, 2.00, 10) data = plot_conductance(sys, Es, 0,1) #case 2 sys, leads = create_system(Norbs=1 ) sys = sys.finalized() kwant.plot(sys,fig_size=(20,20)) Es = np.linspace(0.0001, 2.00, 10) data = plot_conductance(sys, Es, 0,1) if __name__ == "__main__": main() On Mon, Apr 12, 2021 at 11:38 PM Christoph Groth <christoph.groth@cea.fr> wrote:
Henry, Adel, thanks a lot for reporting this problem and finding a workaround.
Setting norbs may be required in some cases, but computation results of Kwant should not depend on whether norbs is set or not! This could be a rather nasty and serious bug, I would like to look closely into it.
Could you show me the exact workaround that you have found? That will hopefully help me to track down this problem.
Cheers Christoph
-- Abbout Adel