A long time ago Maciej Zwierzycki wrote: > 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. ;)
Your problem is not related to a bug, but rather to subtleties of how sites work in Kwant. I will try to explain everything in detail. Perhaps we can come up with ideas how Kwant could be improved to avoid such problems? A site in Kwant always belongs to a site family, i.e. an instance of a class derived from ‘SiteFamiliy’. In fact, a site is just a tuple consisting of a tag and a family to which some custom methods have been added. A site family specify all the properties of the sites that belong to it. A very simple example of a custom site family is shown in the Kwant paper. In practice, almost all site families are instances of the ‘kwant.lattice.Monatomic’ class, which is an abstraction for a Bravais lattice with a single atom basis. It would be possible to create a site family that represents a polyatomic Bravais lattice, but then the tags of such a site family would be a mix of an index for the basis and indices for the actual Bravais lattice. E.g. creating two sites at position (2, 3) that belongs to, respectively, sublattices 0 and 1 would look like this: lat(0, 2, 3) lat(1, 2, 3) Whereas we felt the following to be nicer: lat.a(2, 3) lat.b(2, 3) Note that in the second snippet, ‘lat’ has been created by calling ‘kwant.lattice.honeycomb’ and is an instance of ‘kwant.lattice.Polyatomic’ to which two instance attributes ‘a’ and ‘b’ were added like this lat.a, lat.b = lat.sublattices The object ‘lat’ is *not* a site family. It’s the sublattices that are the site families. ---------------- Two important properties of site families are that they are (effectively) immutable and can be pickled. • Immutability is important so that sites can be used as keys in dictionaries. (This is for example used internally by the builder module, but it’s also handy for users.) • Pickleability is important so that sites and objects containing them (like systems) can be saved to the disk and sent over a network. This comes very handy for parallel computation for example. These two properties have the consequence that two site families that were created with the same parameters are equal. It works just like numbers, or strings, or other immutable objects. >>> lat_a = kwant.lattice.Monatomic([[1, 0], [0, 1]]) >>> lat_b = kwant.lattice.Monatomic([[1, 0], [0, 1]]) >>> lat_a == lat_b True (If one needs multiple separate lattices with the same parameters it’s possible to distinguish by using the ‘name’ parameter.) > 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. ;) If you were to do this with two different monatomic lattices (e.g. a square lattices), you would get an error that the right lead is not stopped by the central region. But in your case the situation is more complicated: Each call to ‘kwant.lattice.general’ returns an instance of ‘Polyatomic’ which is just a collection of monatomic lattices that are the actual site families. In your case, since the first vector of the basis is the same for ‘lat’ and ‘lat2’, but the second differs, we have: >>> lat.sublattices[0] == lat2.sublattices[0] True >>> lat.sublattices[1] == lat2.sublattices[1] False One can see this in the plots where the color of a site depends on the site family. Sublattice ‘lat.sublattices[0]’ (which is equal to ‘lat2.sublattices[0]’) is blue, ‘lat.sublattices[1]’ is orange, and ‘lat2.sublattices[1]’ is green. > 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. The naked eye is misguided because there are two distinct site families that occupy the same real space positions, namely the “orange” and the “green” site family. Kwant is happily plotting several sites on top of each other... The function ‘attach_lead’ is only concerned about sites that belong to the families that are present in the lead. (The real space position of a site does not matter here by design.) So, when attaching the lead from the right, it is only concerned about blue and green sites in the central region. Thus, when filling up the scattering region such that it neatly matches the lead unit cell, ‘attach_lead’ happily adds green sites over already existing orange sites and connects them to the surrounding blue sites. The pre-existing orange sites remain dangling with two hoppings only. I attach a plot created by a slightly modified version of your script that is also attached. Observe that in the plot some green triangles are plotted over orange squares. > 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? Everything is working as designed. However I am always interested to discuss ways in which the design could be improved! > 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? This should be clear from the above explanation now, I hope. > 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? I think that it was a good decision that we decoupled the real space position of a site from its identity. (In fact, real space position for a site is optional in Kwant.) What is debatable is whether it was a wise choice to make sublattices the site families, especially since the user does not create them directly when making a polyatomic lattice. This could be actually changed, I guess even in a practically backwards-compatible way. Any ideas, observations, questions? Cheers Christoph
Figure_2.pdf
Description: Adobe PDF document
import kwant
import numpy as np
from numpy import sin, cos,tan,sqrt
from matplotlib import pyplot
def site_symbol(site):
if site.family in lat.sublattices:
return 'S'
else:
return ('P', 3, 0)
def make_system(a, t, W, L1):
global lat, lat2
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, site_symbol=site_symbol, show=False)
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)
return syst
syst = make_system(a=1.42,t=3.0,W=10.5,L1=20)
kwant.plot(syst, site_symbol=site_symbol)
syst=syst.finalized()
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()
smime.p7s
Description: S/MIME cryptographic signature
