Surprising (?) behaviour of attach leads
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()
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
Christoph Groth wrote:
A long time ago Maciej Zwierzycki wrote:
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?
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.
Here is an idea about how problems similar to yours might be avoided in the future. It seems to me that turning polyatomic lattices into site arrays may be neither necessary nor desirable. Things would have been less confusing for you had you assigned different names to ‘lat’ and ‘lat2’ using the ‘name’ option of ‘kwant.lattice.general’. Then you would have immediately got the error message “Builder does not interrupt the lead, this lead cannot be attached.” However, that ‘name’ option is not widely used and probably often not known. We could make assigning a name to a lattice a requirement. But this would needlessly complicate most simple scripts that happily use a single lattice. So what we could do instead is warn whenever a lattice object is created explicitly that is identical to an already existing one, but unnamed. (Technically, a registry of existing lattices that does not prevent them from being garbage collected can be realized using the weakref module of the standard library.) Christoph
Dear Christoph, Thanks a lot for a very thorough and clear explanation. I think I even understood most of it. :) I do not think there is anything wrong with the current treatment of sublattices, once it is properly understood. It probably boils down to personal taste. Just like thinking in terms of lattice with polyatomic basis vs a collection of monoatomic lattices. Re: "So what we could do instead is warn whenever a lattice object is created explicitly that is identical to an already existing one, but unnamed." Every little helps so why not. On the other hand I suspect the message itself will not be sufficiently clear for someone without the good understanding of how Kwant is representing a physical system/model ... cough ... a Fortran programmer like yours truly. A useful thing would be to actually check for site overlap during attachment of the lead but this might violate your philosophy of decoupling site from its position. I do see where you are coming from. The TB Hamiltonian is indeed well described by the graph structure -- set of nodes plus connections -- with positions important only in the setup phase. On the other hand lead attachment is the setup phase... I guess it all depends on what would be the wording of the warning. A simple statement of there being two identical objects would be sufficient for an experienced user, who might not need it anyway, but useless for a newbie. I have to admit I struggle a bit with the sophistication and the abstraction level on which Kwant operate. But then I never properly learned the whole OO programming paradigm etc... To summarize a warning would be nice, but perhaps it could point to specific problems, at least aas example? In any case thanks again for your answer. It was certainly helpful for me. Regards MZ
participants (3)
-
Christoph Groth -
Maciej Zwierzycki -
maciej.zwierzycki＠ifmpan.poznan.pl