Dear Kwant Authors and Users,
Greetings and thank you for a wonderful tool.
While trying to force Kwant to accept the unit cell of the lead in a particular, connected form, I've encountered somewhat puzzling
situation illustrated in the attached piece of code. The described procedure does not make sense physically, but its purpose is
to illustrate the behaviour which may possibly (?) lead to problems, at least for newbies like me. ;)
Let's say I define two honeycomb lattices, differing in the choice of basis, but entirely equivalent (i.e. they overlap)
prime=np.array([[sqrt(3)/2,sqrt(3)],[1.5,0.0]])*a
basis_at=np.array([[0.0,0.0],[-0.5,0.5]])*a
lat = kwant.lattice.general([prime[:,0],prime[:,1]],[basis_at[:,0],basis_at[:,1]])
and
basis2=np.array([[0.0,sqrt(3)/2],[-0.5,-1.0]])*a
lat2 = kwant.lattice.general([prime[:,0],prime[:,1]],[basis2[:,0],basis2[:,1]])
I'm using numpy so that vector algebra works. The rectangular scattering region and the left lead (same width)
are populated with "lat". Now I'm adding a few atoms from lat2 to the right edge of the scattering region,
e.g. a single hexagon
for i in range (0,1): # changing the range we can generate also the full layer of lat2 sites
syst[lat2.sublattices[0](0+2*i,9-i)]=0.0
syst[lat2.sublattices[1](0+2*i,8-i)]=0.0
syst[lat2.sublattices[0](1+2*i,8-i)]=0.0
syst[lat2.sublattices[1](1+2*i,8-i)]=0.0
and then attach the lead, also populated with "lat2". I know, that's not how it's described in tutorial. ;)
The filling algorithm generates the system which to naked eye looks like the perfect nanowire.
The transmission however is not perfect i.e. T < N (number of modes). If the full layer of lat2 sites is added to the right edge,
than everything is as it should be, i.e. perfect transmission.
I'm not sure if it's a bug - I'm doing things not in the manual - but I find it strange. A perfect system, however generated, should
lead to perfect transmission. What gives?
Also I noticed that neighbours() method called for one of the lattices tends to pick up the sites from the other lattice,
if they belong to he sublattice generated from (0,-0.5a), i.e. from the common basis. Try commenting out one of the lines 53-54.
Is it intentional or not?
It's quite likely that the described situation never arises in real-world situations especially if one does the right thing
and pads the interface with the full principal layer before attaching the leads. I've stumbled upon it only thanks to my laziness.
But it might also illustrate some fragility of the algorithm, possibly worth eliminating.
In any case I'd be grateful for some illumination. Why isn't the "perfect" wire behaving as it should?
Where is the imperfection if it's not visible in the structure?
With Kind Regards
Maciej Zwierzycki
import kwant
import numpy as np
from numpy import sin, cos,tan,sqrt
from matplotlib import pyplot
def make_system(a, t, W, L1):
prime=np.array([[sqrt(3)/2,sqrt(3)],[1.5,0.0]])*a
basis_at=np.array([[0.0,0.0],[-0.5,0.5]])*a
lat = kwant.lattice.general([prime[:,0],prime[:,1]],[basis_at[:,0],basis_at[:,1]])
syst = kwant.Builder()
#### Define the scattering region. ####
def scatt_shape(pos):
(x, y) = pos
return (-L1 < x < L1) and ( -W <y < W)
syst[lat.shape(scatt_shape,basis_at[:,0])]=0
syst[lat.neighbors()] = t
#### Left lead. ####
def lead_shape(pos):
(x, y) = pos
return (-W < y < W )
lsym=kwant.TranslationalSymmetry(lat.vec((0,-1)))
lsym.add_site_family(lat.sublattices[0], other_vectors=[(-2,1)])
lsym.add_site_family(lat.sublattices[1], other_vectors=[(-2,1)])
llead = kwant.Builder(lsym)
llead[lat.shape(lead_shape,basis_at[:,0])] = 0
llead[lat.neighbors()] = t
#### Right lead. ####
### Define a hexagonal lattice with different - but equivalent! - basis
basis2=np.array([[0.0,sqrt(3)/2],[-0.5,-1.0]])*a
lat2 = kwant.lattice.general([prime[:,0],prime[:,1]],[basis2[:,0],basis2[:,1]])
#### Extra hexagon of lat2 on the right edge of the scattering region:
#### -- range(0,1) - single hexagon
##### -- range(-2,3) - full width
for i in range (0,1):
syst[lat2.sublattices[0](0+2*i,9-i)]=0.0
syst[lat2.sublattices[1](0+2*i,8-i)]=0.0
syst[lat2.sublattices[0](1+2*i,8-i)]=0.0
syst[lat2.sublattices[1](1+2*i,8-i)]=0.0
syst[lat2.neighbors()]=t
syst[lat.neighbors()]=t
kwant.plot(syst)
lsym=kwant.TranslationalSymmetry(lat2.vec((0,1)))
lsym.add_site_family(lat2.sublattices[0], other_vectors=[(2,-1)])
lsym.add_site_family(lat2.sublattices[1], other_vectors=[(2,-1)])
rlead= kwant.Builder(lsym)
rlead[lat2.shape(lead_shape,basis_at[:,0])] = 0
rlead[lat2.neighbors()] =t
syst.attach_lead(llead)
syst.attach_lead(rlead)
syst=syst.finalized()
return syst
syst = make_system(a=1.42,t=3.0,W=10.5,L1=20)
kwant.plot(syst)
data1 = []
data2 = []
params = dict(a=1.42, t=-3.0)
energies=[0.2*i+0.01 for i in range(20)]
for energy in energies:
print('En=',energy)
smatrix = kwant.smatrix(syst, energy,params=params)
data1.append(smatrix.transmission(1, 0))
data2.append(smatrix.num_propagating(0))
pyplot.figure()
pyplot.plot(energies, data1, 'r--', label='T')
pyplot.plot(energies, data2,'r>', label='N')
legend = pyplot.legend(loc='upper left')
pyplot.xlabel("energy")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()