Hello, I am having trouble figuring out what the cause of the change in conductance is when the only thing I do is I add the "norbs=1" command to the definition of the electron and hole lattices (is this not default?). Below is the code, I have added in "norbs = 1" in lines 70 and 71, if you take them out, you get the proper result. (sorry of the long code). 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_length, np_transition_length, p_length, width, lattice_spacing, onsite_potential_left_lead, onsite_potential_superconducting_lead, onsite_potential_n, onsite_potential_p, hopping_parameter, Delta, boundary_hopping): saved_args = locals() padding = 0.5*lattice_spacing*tan_30 def calc_total_length(n_length, np_transition_length, p_length, lattice_spacing): total_length = n_length + np_transition_length + p_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: n_length = n_length - diff total_length = new_length return total_length, n_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 potential_e(site): (x,y) = site.pos return onsite_potential_n + (0.5*(1+tanh((x-n_length- 0.5*np_transition_length)/np_transition_length)))*(onsite_potential_p- onsite_potential_n) def potential_h(site): (x,y) = site.pos return (-1)*(onsite_potential_n + (0.5*(1+tanh((x-n_length- 0.5*np_transition_length)/np_transition_length)))*(onsite_potential_p- onsite_potential_n)) 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)) total_length, n_length = calc_total_length(n_length, np_transition_length,p_length, lattice_spacing) width, N = calc_width(width, lattice_spacing, padding) ### Defining lattices ### graphene_e = kwant.lattice.honeycomb(a=lattice_spacing,name='e', norbs=1) graphene_h = kwant.lattice.honeycomb(a=lattice_spacing,name='h', norbs=1) a, b = graphene_e.sublattices c, d = graphene_h.sublattices sys = kwant.Builder() sys[graphene_e.shape(rectangle, (0.5*lattice_spacing, 0))] = potential_e sys[graphene_h.shape(rectangle, (0.5*lattice_spacing, 0))] = potential_h sys[graphene_e.neighbors()] = -hopping_parameter sys[graphene_h.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 if str(site.family) == "<Monatomic lattice h1>": sys[d(x,y),c(int(x+tag_site_calc(N)),N)] = boundary_hopping # kwant.plot(sys) kwant.plot(sys, site_color=potential_e, site_lw=0,hop_lw=0,site_size=0.2, lead_site_lw=0, colorbar=True, cmap='jet') ### Defning leads ### # Symmetry of the left and right leads sym_l = kwant.TranslationalSymmetry((-1, 0)) sym_r = kwant.TranslationalSymmetry((1, 0)) # Redefining the shape of the unit cell so that we don't add new sites sym_l.add_site_family(graphene_e.sublattices[0], other_vectors=[(-1, 2)]) sym_l.add_site_family(graphene_e.sublattices[1], other_vectors=[(-1, 2)]) sym_r.add_site_family(graphene_e.sublattices[0], other_vectors=[(-1, 2)]) sym_r.add_site_family(graphene_e.sublattices[1], other_vectors=[(-1, 2)]) sym_l.add_site_family(graphene_h.sublattices[0], other_vectors=[(-1, 2)]) sym_l.add_site_family(graphene_h.sublattices[1], other_vectors=[(-1, 2)]) sym_r.add_site_family(graphene_h.sublattices[0], other_vectors=[(-1, 2)]) sym_r.add_site_family(graphene_h.sublattices[1], other_vectors=[(-1, 2)]) lead0 = kwant.Builder(sym_l) # Left electron lead0[graphene_e.shape(lead_shape, (1,1))] = -onsite_potential_left_lead lead0[graphene_e.neighbors()] = -hopping_parameter for site in lead0.sites(): (x,y) = site.tag if site.pos[1]<0: lead0[b(x,y),a(int(x+tag_site_calc(N)),N)] = -boundary_hopping lead1 = kwant.Builder(sym_l) # Left hole lead1[graphene_h.shape(lead_shape, (1,1))] = onsite_potential_left_lead lead1[graphene_h.neighbors()] = hopping_parameter for site in lead1.sites(): (x,y) = site.tag if site.pos[1]<0: lead1[d(x,y),c(int(x+tag_site_calc(N)),N)] = boundary_hopping lead3 = kwant.Builder(sym_r) # Right electron lead3[graphene_e.shape(lead_shape, (1,1))] = -onsite_potential_superconducting_lead lead3[graphene_e.neighbors()] = -hopping_parameter for site in lead3.sites(): (x,y) = site.tag if site.pos[1]<0: lead3[b(x,y),a(int(x+tag_site_calc(N)),N)] = -boundary_hopping lead4 = kwant.Builder(sym_r) # Right hole lead4[graphene_h.shape(lead_shape, (1,1))] = onsite_potential_superconducting_lead lead4[graphene_h.neighbors()] = hopping_parameter for site in lead4.sites(): (x,y) = site.tag if site.pos[1]<0: lead4[d(x,y),c(int(x+tag_site_calc(N)),N)] = boundary_hopping 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)] = Delta if str(site.family) == "<Monatomic lattice e1>": lead2[d(x,y),b(x,y)] = Delta # 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], width, N, saved_args def plot_conductance(syst, energies, lead_in, lead_out, saved_args, width, geo_prop, pot_prop): """ Return conductance as a function of energy from lead_in to lead_out. """ lead_dict = { 0: "Electron Lead", 1: "Hole Lead", 2: "Superconducting Lead" } data = [] step_counter = 1 energy_size = energies.size elapsed_time_per_step_average = 0.0 for energy in energies: start = time.time() print('Step: ', step_counter, '/', energy_size) smatrix = kwant.smatrix(syst, energy) data.append(smatrix.transmission(lead_out, lead_in)) end = time.time() elapsed_time_per_step = int(end - start) print(' Time per step:', elapsed_time_per_step, ' seconds') if elapsed_time_per_step_average > 0: elapsed_time_per_step_average = (elapsed_time_per_step_average + elapsed_time_per_step)/2 else: elapsed_time_per_step_average = elapsed_time_per_step time_left = int(elapsed_time_per_step_average*(energy_size-step_counter)/60) print(' Time left: ', time_left, ' minutes') step_counter = step_counter + 1 pyplot.figure() pyplot.plot(energies, data) str1 = str(lead_dict[lead_in]) str2 = str(lead_dict[lead_out]) title_str = "From " + str1 + " to " + str2 pyplot.title(title_str) pyplot.xlabel("Energy [t]") pyplot.ylabel("Conductance [e^2/h]") ### All variables used ### # comment2_txt = "N Length = " + str(saved_args['n_length']) + ", NP Length = " + str(saved_args['np_transition_length']) + ", P Length = " + str(saved_args['p_length']) + ", Width = " + str(0.01*int(100*width)) + ", Potential Left Lead = " + str(saved_args['onsite_potential_left_lead']) + ", Potential Superconducting Lead = " + str(saved_args['onsite_potential_superconducting_lead']) + ", Potential N = " + str(saved_args['onsite_potential_n']) + ", Potential P = " + str(saved_args['onsite_potential_p']) + ", Hopping Parameter = " + str(saved_args['hopping_parameter']) + ", Delta = " + str(saved_args['Delta']) + ", Resolution = " + str(energy_size) ### Main variables: pn_transition_length, width, length_of_P ### comment2_txt = "NP Transition Length = " + str(geo_prop[1]) + ", P Length = " + str(geo_prop[2]) + ", Width = " + str(geo_prop[0]) comment1_txt = "Left Lead = " + str(pot_prop[0]) + ", N-Region: " + str(pot_prop[1]) + ", P-Region: "+str(pot_prop[2]) +", SC Lead: " + str(pot_prop[3])+", Delta: " +str(pot_prop[4]) comment0_txt = "Resoltuion: " + str(energy_size) fig_txt = tw.fill(tw.dedent(comment1_txt.rstrip() ), width=80) pyplot.figtext(0.5, -0.12, fig_txt, horizontalalignment='center', fontsize=8, multialignment='left', bbox=dict(boxstyle="round", facecolor='#D8D8D8', ec="0.5", pad=0.5, alpha=1), fontweight='bold') fig_txt = tw.fill(tw.dedent(comment2_txt.rstrip() ), width=80) pyplot.figtext(0.5, -0.05, fig_txt, horizontalalignment='center', fontsize=8, multialignment='left', bbox=dict(boxstyle="round", facecolor='#D8D8D8', ec="0.5", pad=0.5, alpha=1), fontweight='bold') fig_txt = tw.fill(tw.dedent(comment0_txt.rstrip() ), width=80) pyplot.figtext(0.81, 0.84, fig_txt, horizontalalignment='center', fontsize=6, multialignment='left', bbox=dict(boxstyle="round", facecolor='#D8D8D8', ec="0.5", pad=0.5, alpha=1), fontweight='bold') # pyplot.ylim(0,3) pyplot.xlim(0.0,2.0) pyplot.show() return data def calc_fermi_wavelength(onsite_potential): return 2*pi*(1/sqrt(3+onsite_potential)) def main(): """ GEOMETRIC PROPERTIES - IN UNITS OF FERMI WAVELENGTH """ geo_prop = np.array([8, # Width 2, # PN-Junction transition length 2]) # P-Region length """ ENERGIES - IN UNITS OF t """ pot_prop = np.array([-1.0001, # Left lead 0.3, # N-Region 0.3, # P-Region 0.3, # Superconducting lead 1.0]) # Delta # Converting property arrays from units of Fermi wavelength to #units of lattice sites. fermi_wavelength = calc_fermi_wavelength(pot_prop[2]) geo_prop2 = fermi_wavelength*geo_prop k_phase_kick = 0.0 # Phase-kick due to boundary condition ## create_system using above defined property arrays sys, leads, width, N, saved_args = create_system(n_length = ceil(2*fermi_wavelength), np_transition_length = geo_prop[1], p_length = ceil(geo_prop2[2]) # GEOMETRIC PROPERTIES , width = geo_prop2[0], lattice_spacing = 1.0, onsite_potential_left_lead = pot_prop[0], onsite_potential_superconducting_lead = pot_prop[3], # POTENTIALS onsite_potential_n = pot_prop[1], onsite_potential_p = pot_prop[2], hopping_parameter = 1.0, Delta = pot_prop[4], # HOPPING PARAMETERS boundary_hopping = 1.0*exp(-1j*k_phase_kick)) sys = sys.finalized() ### Plot Conductance ### Es = np.linspace(0.0001, 2.00, 20) data = plot_conductance(sys, Es, 0,1, saved_args, width, geo_prop, pot_prop) return #d_energies if __name__ == "__main__": main()