Adding "norbs = 1" changes conductance calculation
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()
Dear Henry, Indeed, when you put norbs=1, you get the wrong result. I want just to highlight that "norbs=1" is not the default value. You can check this in the source code (the default value is "norbs=None"). In my understanding, for "norbs=1", the program is expecting from you to put a (1,1) array as a potential rather than a scalar. This should be the source of your error. @Chridtoph_Groth, I don't know if this is something that needs to be highlighted in the documentation. I hope this helps. Adel On Fri, Apr 9, 2021 at 3:35 PM Henry Axt <henry.axt@gmail.com> wrote:
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()
-- Abbout Adel
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
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
Dear Henry, The problem you experience comes from the fact, that if you specify `norbs` name of the lattice changes from '<Monatomic lattice e0>' to '<Monatomic lattice e0 with 1 orbitals>' while you conditionally insert hoppings, comparing site family string representation with the first variant. The best approach to fix this is to compare site families directly and not with their string representation, this will be faster and more reliable. With best regards, Viacheslav Ostroukh P.S. @Abbout Adel: special thanks for cleaning up the code, helped a lot! On 16/04/2021 02:52, Abbout Adel wrote:
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 <mailto: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
-- With best regards, Slava Ostroukh
Dear all, I think that selecting site families by name is the most clear: `if site.family.name == 'e': ...`. On the other hand, defining `norbs=None` vs. `norbs=1` makes that the onsite values are broadcasted to a tinyarray of shape `(1,)`. This might lead to uncommon broadcasting issues, and data type conversion such as `int` to `float`. Best regards, Pablo
Hi, Site families are comparable and even hashable, and are quite efficient at this, so there’s actually nothing wrong with comparing them directly. So instead of Henry’s (that involves formatting a new string each time) if str(site.family) == "<Monatomic lattice e0>": ... one could use the clearer, more robust, and more efficient if site.family == a: ... The above uses the fact that the sublattices of graphene_e were given names: a, b = graphene_e.sublattices Cheers Christoph
participants (5)
-
Abbout Adel
-
Christoph Groth
-
Henry Axt
-
pablo.perez.piskunow@gmail.com
-
Viacheslav Ostroukh