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()