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.karbaschi@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