Hello All,

I am trying to simulate black phosphorous quantum dots of circular shape, but with their ends (that will be attached to the leads) will be broad. The shape would be something like a circle with two thick lines on either side of the circle. I have the following code for it, but the problem I face is that the two thick lines are disconnected from the circle. I am not able to understand the reason behind it and rectify it. If any of you could help me, it would be of great value to me.

------ CODE -------

import kwant
from matplotlib import pyplot
import numpy as np

h = 2.707 * np.cos(13.69*np.pi/180)
dx = 4.43
dy = 3.27
shift_x = 2.164 * np.cos(0.5*98.15*np.pi/180)
shift_y = 2.164 * np.sin(0.5*98.15*np.pi/180)

r_const = 20

lat0 = kwant.lattice.general([(dx, 0, 0), (0, dy, 0)], [(0, 0, 0), (shift_x, shift_y, 0)])
lat1 = kwant.lattice.general([(dx, 0, 0), (0, dy, 0)], [(0.5*dx, 0.5*dy, h), (0.5*dx + shift_x, 0.5*dy + shift_y, h)])

def cuboid_shape(pos):
x, y, z = pos
if -1.2*r_const <= x < -r_const:
return abs(y) < r_const
if -r_const <= x <= r_const:
return x**2 + y**2 <= r_const**2
if r_const < x <= 1.2*r_const:
return abs(y) < r_const

def potential(site):
x, y, z = site.pos
return (-tanh(5*(x-r_const)) - tanh(5*(x+r_const))) * pot

sys = kwant.Builder()
sys[lat0.shape(cuboid_shape, (0, 0, 0))] = potential
sys[lat0.neighbors(1)] = 1
sys[lat0.neighbors(2)] = 1
sys[lat1.shape(cuboid_shape, (0.5*dx, 0.5*dy, h))] = potential
sys[lat1.neighbors(1)] = 1
sys[lat1.neighbors(2)] = 1

lat0_0 = [ ]
lat0_1 = [ ]
lat1_0 = [ ]
lat1_1 = [ ]

for site in kwant.Builder.sites(sys):
if site.family == lat0.sublattices[0]:
lat0_0.append(site[1])
elif site.family == lat0.sublattices[1]:
lat0_1.append(site[1])
elif site.family == lat1.sublattices[0]:
lat1_0.append(site[1])
elif site.family == lat1.sublattices[1]:
lat1_1.append(site[1])

for w in lat0_0:
for v in lat0_1:
for k in lat1_0:
for l in lat1_1:
if v == k:
sys[lat0.sublattices[1](v[0],v[1]), lat1.sublattices[0](k[0],k[1])] = 2
if w == l:
if (w[0]+1,w[1]+1) in lat0_0:
sys[lat0.sublattices[0](w[0]+1,w[1]+1), lat1.sublattices[1](l[0],l[1])] = 2
if w == l:
if (l[0]-1,l[1]-1) in lat1_1:
sys[lat0.sublattices[0](w[0],w[1]), lat1.sublattices[1](l[0]-1,l[1]-1)] = 2

kwant.plot(sys, site_size=0.2, hop_lw=0.1)

Any help would be greatly appreciated.

Regards,
Shivang Agarwal

--
Shivang Agarwal