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) On Fri, Aug 16, 2019 at 10:13 PM Hadi Zahir <zahir.sci@gmail.com> wrote:
Thanks, Anton!
Creating a fake system is a good idea, but it seems that it is still not a finalized system as it returns the following error:
"TypeError: The system needs to be finalized."
Here is the piece of code I use:
syst_c=[] for item in finilized_sys.graph: if finalized_sys.sites[item[0]].family in [a1, b1] : syst_c.append syst_slice = SimpleNamespace(sites=syst_c) densities, box = kwant.plotter.interpolate_density(syst=syst_slice, density=densities)
Best, Hadi
On Friday, August 16, 2019, 9:16:29 AM GMT+1, Anton Akhmerov < anton.akhmerov+kd@gmail.com> wrote:
Hi Hadi,
The intended way would be to extend the density vector with np.nan, so that it has a correct length. However since you mention performance, I'm not sure if this would achieve what you are looking for. Inspecting the source of interpolate_density [1], we see that the only way the system is used is to get
sites = np.array([s.pos for s in syst.sites])
Therefore we can mock the necessary properties of the system by using the following code:
from types import SimpleNamespace # make a fake system that only contains the correct sites and isn't actually a system syst_slice = SimpleNamespace(sites=(site for site in syst.sites if where(site)))
It is certainly a hack, so a more reliable solution would be to take over the bits of the source of interpolate_density and the low level function that it uses [2] and re-purpose for your needs.
[1]: https://gitlab.kwant-project.org/kwant/kwant/blob/v1.4.1/kwant/plotter.py#L1... [2]: https://gitlab.kwant-project.org/kwant/kwant/blob/v1.4.1/kwant/plotter.py#L1...
On Fri, 16 Aug 2019 at 10:03, Hadi Zahir <zahir.sci@gmail.com> wrote:
Dear Kwant developers,
I have a finalized closed system ''sys_f'' for which I computed the
local charge densities "ch_densities", which is basically a vector of length=len(sys_f.sites).
For complexity reason, I'd like to feed a slice of sys_f to the
interpolate_density() as:
kwant.plotter.interpolate_density(syst=sys_p, density=ch_densities_p)
where sys_p is a slice of sys_f and ch_densities_p is the corresponding
density vector.
Could you please let me know how to do it?
Best regards, Hadi
-- Abbout Adel