Dear Zahir,
You can for example make a copy of your initial system and then delete all the region you do not want to keep.
You can do that as follows:
sys_p=kwant.update(syst)
It is important that you get also the indices of the sites you kept in order to extract the density of the remaining region from the initial density array:
Density_p=Density[indices] #here indices is a list of the indices of the sites you kept.
It is also important to add the leads in your original system only after you do the copy of your system and delete the unwanted sites.
Now, depending on what you need, you can feed your function with sys_p or sys_p.finalized().
Check the example below for more details.
I hope this helps.
Regards,
Adel
####################################################################
from math import pi, sqrt, tanh
import numpy as np
import kwant
# For computing eigenvalues
import scipy.sparse.linalg as sla
# For plotting
from matplotlib import pyplot
# Define the graphene lattice
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
[(0, 0), (0, 1 / sqrt(3))])
a, b = graphene.sublattices
def make_system(r=10, w=2.0, pot=0.1):
#### Define the scattering region. ####
# circular scattering region
def circle(pos):
x, y = pos
return x ** 2 + y ** 2 < r ** 2
syst = kwant.Builder()
# w: width and pot: potential maximum of the p-n junction
def potential(site):
(x, y) = site.pos
d = y * cos_30 + x * sin_30
return pot * tanh(d / w)
syst[graphene.shape(circle, (0, 0))] = potential
# specify the hoppings of the graphene lattice in the
# format expected by builder.HoppingKind
hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
# Modify the scattering region
del syst[a(0, 0)]
syst[a(-2, 1), b(2, 2)] = -1
#### Define the leads. ####
# left lead
sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
def lead0_shape(pos):
x, y = pos
return (-0.4 * r < y < 0.4 * r)
lead0 = kwant.Builder(sym0)
lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
# The second lead, going to the top right
sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))
def lead1_shape(pos):
v = pos[1] * sin_30 - pos[0] * cos_30
return (-0.4 * r < v < 0.4 * r)
lead1 = kwant.Builder(sym1)
lead1[graphene.shape(lead1_shape, (0, 0))] = pot
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
return syst, [lead0, lead1]
pot = 0.1
syst, leads = make_system(pot=pot)
# To highlight the two sublattices of graphene, we plot one with
# a filled, and the other one with an open circle:
def family_colors(site):
return 0 if site.family == a else 1
# Plot the closed system without leads.
kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False)
# In[ ]:
Sites=list(syst.finalized().sites)
# In[ ]:
def region(pos):
x,y=pos
return .2*y**2+abs(3*x)<15
def color(site):
if region(site.pos): return 'g'
else: return 'k'
#this function makes a copy of the system
# and returns the indices of the sites in the region you want to keep
def Fun(syst,point):
sys_p=kwant.Builder()
sys_p.update(syst) #make a copy of the system
sites_to_keep=[site for site in Sites if region(site.pos)]
sites_to_del=[site for site in Sites if (site not in sites_to_keep)]
indices=[syst.finalized().id_by_site[site] for site in sites_to_keep]
del sys_p[sites_to_del]
return sys_p,indices
kwant.plot(syst,site_color=color)
# In[ ]:
sys_p,indices=Fun(syst,(0,1))# important to do this before attaching the leads
for lead in leads: syst.attach_lead(lead)
wave=kwant.wave_function(syst.finalized(),energy=0.3)
wf=wave(0)[0]
density=abs(wf)**2
# In[ ]:
Sites_p=list(sys_p.finalized().sites)
density_p=density[indices]
def size(i):
return 5.5*density_p[i]/density.max()
kwant.plotter.map(sys_p.finalized(),density_p)
kwant.plot(sys_p.finalized(),site_size=size,site_color=(0, 0, 1, 0.6))
# In[ ]:
kwant.plot(syst,site_color=color)