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