Dear Hossein,
You can do that easily with "del":
import random
Sites=syst.finalized().sites
To_delete=random.sample(Sites, 23) # choose 23 sites to delete
del syst[To_delete]
kwant.plot(syst)
I hope this helps,
Adel
On Tue, Jun 20, 2023 at 3:47 PM Hossein Karbaschi <h.karbas...@gmail.com>
wrote:
> Hi everybody,
>
> In the following code I want to delete few random sites from the
> scattering region.
> Please help me in this case.
> Best regards,
> Hossein
>
> import kwant
> from math import sqrt
> import matplotlib.pyplot as plt
> import tinyarray
> import numpy as np
> import math
>
>
> import matplotlib
>
>
> #Constants.........................................................
>
> #Parameter.........................................................
> a0=2.617*sqrt(3);
> d0=2.617;
> dz0=1.59047;
> D0=sqrt((d0)**2+(dz0)**2);
>
> #st=0;
> st=-0.0;
> dz=1.59047+st;
> d=-0.3*dz+3.09;
> a1=-0.3*sqrt(3)*dz+5.35;
> D=sqrt((d)**2+(dz)**2);
>
> scale=(D/D0)**(-8);
> qo=[];
> N=dz/D;
> delta1=0j;
> delta2=0;
> R=(2*a1)**2;
> #on-site energy...................................................
> Es=-10.906;
> Ep=-0.486;
>
> Ep1=-0.486;
> Ep2=-0.486;
>
> #nearest neighbor................................................
> sssigma=-0.608*scale;
> spsigma=1.32*scale;
> ppsigma=1.854*scale;
> pppay=-0.6*scale;
>
> #next nearest neighbor.........................................
> ppsigma2=0.156*scale;
> pppay2=0;
> sssigma2=0;
> spsigma2=0;
>
> #onsite potentials ...........................................
> pot_a=10;
> pot_b=0;
> #nearest
> neighbors.............................................................
> l1=(math.cos(math.pi/6))*(d/D);
> m1=(math.sin(math.pi/6))*(d/D);
> N1=dz/D;
>
> #second
> neighbors..............................................................
> l2=(math.cos(math.pi/3))
> m2=(math.sin(math.pi/3))
>
> latt = kwant.lattice.general([(a1,0),(a1*0.5,a1*math.sqrt(3)/2)],
> [(a1/2,-d/2),(a1/2,d/2)])
>
>
> a,b = latt.sublattices
> syst= kwant.Builder()
>
> # Onsite matrice energt
> Onsite_a=tinyarray.array([[Es,0,0,0],[0,Ep,delta1,0],[0,-delta1,Ep,0],
> [0,0,0,Ep1]])
>
> Onsite_b=tinyarray.array([[Es,0,0,0],[0,Ep,delta1,0],[0,-delta1,Ep,0],
> [0,0,0,Ep2]])
>
>
> pot_edge_a=tinyarray.array([[Es+pot_a,0,0,0],[0,Ep+pot_a,delta1,0],[0,-delta1,Ep+pot_a,0],
> [0,0,0,Ep1+pot_a]])
>
>
> pot_edge_b=tinyarray.array([[Es+pot_b,0,0,0],[0,Ep+pot_b,delta1,0],[0,-delta1,Ep+pot_b,0],
> [0,0,0,Ep2+pot_b]])
> #nearest
> neighbors.............................................................
> first_up=tinyarray.array([[sssigma, l1*spsigma,m1*spsigma,N1*spsigma],
> [-l1*spsigma,
> l1**2*ppsigma+(1-l1**2)*pppay,l1*m1*(ppsigma-pppay),l1*N1*(ppsigma-pppay)],
> [-m1*spsigma ,l1*m1*(ppsigma-pppay),
> m1**2*ppsigma+(1-m1**2)*pppay,N1*m1*(ppsigma-pppay)],
>
> [-N1*spsigma,l1*N1*(ppsigma-pppay),N1*m1*(ppsigma-pppay),N1**2*ppsigma+(1-N1**2)*pppay]])
>
> first_left_up=tinyarray.array([[sssigma,-l1*spsigma,m1*spsigma,N1*spsigma],
> [l1*spsigma, l1**2*ppsigma+(1-l1**2)*pppay,-
> l1*m1*(ppsigma-pppay),-l1*N1*(ppsigma-pppay)],
> [-m1*spsigma ,-l1*m1*(ppsigma-pppay),
> m1**2*ppsigma+(1-m1**2)*pppay,N1*m1*(ppsigma-pppay)],
>
> [-N1*spsigma,-l1*N1*(ppsigma-pppay),N1*m1*(ppsigma-pppay),N1**2*ppsigma+(1-N1**2)*pppay]])
>
> first=tinyarray.array([[sssigma,0,-(d/D)*spsigma,N1*spsigma],
> [0,pppay,0,0],
> [(d/D)*spsigma ,0,
> (d/D)**2*ppsigma+(1-(d/D)**2)*pppay,-N1*(d/D)*(ppsigma-pppay)],
>
> [-N1*spsigma,0,-N1*(d/D)*(ppsigma-pppay),N1**2*ppsigma+(1-N1**2)*pppay]])
> #second
> neighbors..............................................................
> second_up=tinyarray.array([[sssigma2,l2*spsigma2,m2*spsigma2,0],
>
> [-l2*spsigma2,l2**2*ppsigma2+(1-l2**2)*pppay2,l2*m2*(ppsigma2-pppay2),0],
> [-m2*spsigma2
> ,l2*m2*(ppsigma2-pppay2),m2**2*ppsigma2+(1-m2**2)*pppay2,0],
> [0,0,0,pppay2]])
>
> second_left_up=tinyarray.array([[sssigma2,-l2*spsigma2,m2*spsigma2,0],
>
> [l2*spsigma2,l2**2*ppsigma2+(1-l2**2)*pppay2,-l2*m2*(ppsigma2-pppay2),0],
> [-m2*spsigma2
> ,-l2*m2*(ppsigma2-pppay2),m2**2*ppsigma2+(1-m2**2)*pppay2,0],
> [0,0,0,pppay2]])
>
> second_right_down=tinyarray.array([[sssigma2,l2*spsigma2,-m2*spsigma2,0],
>
> [-l2*spsigma2,l2**2*ppsigma2+(1-l2**2)*pppay2,-l2*m2*(ppsigma2-pppay2),0],
>
> [m2*spsigma2,-l2*m2*(ppsigma2-pppay2),m2**2*ppsigma2+(1-m2**2)*pppay2,0],
> [0,0,0,pppay2]])
>
> second_down=tinyarray.array([[sssigma2,-l2*spsigma2,-m2*spsigma2,0],
>
> [l2*spsigma2,l2**2*ppsigma2+(1-l2**2)*pppay2,l2*m2*(ppsigma2-pppay2),0],
> [m2*spsigma2
> ,l2*m2*(ppsigma2-pppay2),m2**2*ppsigma2+(1-m2**2)*pppay2,0],
> [0,0,0,pppay2]])
>
> second_right=tinyarray.array([[sssigma2,spsigma2,0,0],
> [-spsigma2,ppsigma2,0,0],[0 ,0,pppay2,0],
> [0,0,0,pppay2]])
>
> second_left=tinyarray.array([[sssigma2,-spsigma2,0,0],
> [spsigma2,ppsigma2,0,0],[0 ,0,pppay2,0],
> [0,0,0,pppay2]])
>
>
> #...................................................................................
> def rectangle(pos):
> x, y = pos
> # z=x**2+y**2
> return -9.2*a1<x<9.2*a1 and -12*d<=y<12*d
>
>
>
> syst[a.shape(rectangle, (1,1))] = Onsite_a
> syst[b.shape(rectangle, (1,1))] = Onsite_b
>
>
> #nearest
> neighbors.............................................................
> syst[[kwant.builder.HoppingKind((0,0),a,b)]] = first
> syst[[kwant.builder.HoppingKind((0,1),a,b)]] = first_up
> syst[[kwant.builder.HoppingKind((-1,1),a,b)]] = first_left_up
>
> #second
> neighbors..............................................................
> syst[[kwant.builder.HoppingKind((0,1),a,a)]]=second_up
> syst[[kwant.builder.HoppingKind((-1,1),a,a)]]=second_left_up
> syst[[kwant.builder.HoppingKind((1,-1),a,a)]]=second_right_down
> syst[[kwant.builder.HoppingKind((0,-1),a,a)]]=second_down
> syst[[kwant.builder.HoppingKind((1,0),a,a)]]=second_right
> syst[[kwant.builder.HoppingKind((-1,0),a,a)]]=second_left
>
> syst[[kwant.builder.HoppingKind((0,1),b,b)]]=second_up
> syst[[kwant.builder.HoppingKind((-1,1),b,b)]]=second_left_up
> syst[[kwant.builder.HoppingKind((1,-1),b,b)]]=second_right_down
> syst[[kwant.builder.HoppingKind((0,-1),b,b)]]=second_down
> syst[[kwant.builder.HoppingKind((1,0),b,b)]]=second_right
> syst[[kwant.builder.HoppingKind((-1,0),b,b)]]=second_left
>
>
>
> sym = kwant.TranslationalSymmetry(latt.vec((1,0)))
>
> sym.add_site_family(latt.sublattices[0], other_vectors=[(-1, 2)])
> sym.add_site_family(latt.sublattices[1], other_vectors=[(-1, 2)])
> lead = kwant.Builder(sym)
>
>
> # ---- random impurity
>
>
>
>
> # --- leads
> def lead_shape(pos):
> x, y = pos
> return -12*d<=y<12*d and -5.2*a1<x<5.2*a1
>
> lead[a.shape(lead_shape, (1,1))] = Onsite_a
> lead[b.shape(lead_shape, (1,1))] = Onsite_b
>
> lead[[kwant.builder.HoppingKind((0,0),a,b)]] =first
> lead[[kwant.builder.HoppingKind((0,1),a,b)]] =first_up
> lead[[kwant.builder.HoppingKind((-1,1),a,b)]] =first_left_up
>
> lead[[kwant.builder.HoppingKind((0,1),a,a)]]=second_up
> lead[[kwant.builder.HoppingKind((-1,1),a,a)]]=second_left_up
> lead[[kwant.builder.HoppingKind((1,-1),a,a)]]=second_right_down
> lead[[kwant.builder.HoppingKind((0,-1),a,a)]]=second_down
> lead[[kwant.builder.HoppingKind((1,0),a,a)]]=second_right
> lead[[kwant.builder.HoppingKind((-1,0),a,a)]]=second_left
>
> lead[[kwant.builder.HoppingKind((0,1),b,b)]]=second_up
> lead[[kwant.builder.HoppingKind((-1,1),b,b)]]=second_left_up
> lead[[kwant.builder.HoppingKind((1,-1),b,b)]]=second_right_down
> lead[[kwant.builder.HoppingKind((0,-1),b,b)]]=second_down
> lead[[kwant.builder.HoppingKind((1,0),b,b)]]=second_right
> lead[[kwant.builder.HoppingKind((-1,0),b,b)]]=second_left
>
> syst.attach_lead(lead,add_cells=0)
>
>
> syst.attach_lead(lead.reversed(),add_cells=0)
> fsys = syst.finalized()
>
>
>
> def family_color(site):
> #print(sys[site])
> if syst[site]==Onsite_a:
> return 'black'
> elif syst[site]==Onsite_b:
> return 'green'
> elif syst[site]==pot_edge_a:
> return 'yellow'
> else:
> return 'blue'
>
> ax=kwant.plot(syst,site_color=family_color);
> ax.savefig('syst.pdf')
>
--
Abbout Adel
