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 <[email protected]> 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
