Dear Sir,
I am a PhD student of Hong Kong University of Science and Technology. I
want to use KWANT to caculate Hall resistance of a Hall bar structure.We
can get the conductance between 6 electrodes, but how to get hall
resistance? Can you give me some help? Thank you very much.
Best Regards,
Zhang Bing

Dear all,
I am using Kwant to study a system with both spins and electron and hole,
so the hopping matrix is 4X4. Now I want to know how to separate the spins,
electron and hole. I study this through the "2.6. Superconductors: orbital
degrees of freedom, conservation laws and symmetries"
smatrix = kwant.smatrix(syst, energy)
data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0]
-smatrix.transmission((0, 0), (0, 0)) +smatrix.transmission((0, 1), (0, 0)))
I do not find more information about this for kwant.
For my system, if I have 3 leads, and the hopping and onsite energy is set
in this order:
(e↑,0,0,0)- spin up electron,first row of the 4X4matrix
(0,e↓,0,0)- spin down electron,second row of the 4X4matrix
(0,0, h↑ ,0)- spin up hole,third row of the 4X4matrix
(0,0, 0, h ↓ )- spin down hole,fourth row of the 4X4matrix
I want to obtain the transmissions from 0 →2.
smatrix = kwant.smatrix(syst, energy)
The transmission from h↑ to e↑ is: smatrix.transmission((2, 1), (0, 2))
The transmission from h ↓ to e↑ is: smatrix.transmission((2, 1), (0, 3))
Is my understanding correct?
Thanks very much in advance!
Hosein Khani

Hi, colleagues!
I'm defining the phosphorene unit cell perfectly. However, when the leads
are attached to the scattering region, the kwant increases the unit cell,
and thus "duplicate" energy levels appear. I'm using the model with 10
interlayer hoppings and 5 interlayer hoppings. Apparently, in the model of
5 interlayer hoppings it doesn't happen! Thank you for your help now.

Hi,
I am trying to calculate the conductance through coupled quantum nanowires.
I have basically considered a square lattice with hopping tx in x direction
and modulating the coupling between the wires through ty i.e. the hopping
in y direction. My system requires periodic boundary condition along x
direction and open boundary condition in y direction. Can you please help
me in how this can be done using kwant to meet the above boundary
conditions. I have mentioned my code below. Since I am new to Kwant, it
would be of much help to me if you help me solving this.
Regards
Deepti Rana
def make_sys1(a=1,W=4,barrierpos=1):
def onsite_normal(site,p):
return( (2 * (p.tx+p.ty) - p.mu) * pauli.s0sz+p.Ez*pauli.szs0))
def onsite_sc(site,p):
return (2 * (p.tx+p.ty ) *
pauli.s0sz+p.Ez*pauli.szs0)+p.delta*pauli.s0sx)
def onsite_barrier(site,p):
return((2*(p.tx+p.ty)+ p.Vbarrier-p.mu)*pauli.s0sz+p.Ez*pauli.szs0))
def hopx(site1, site2, p):
return -p.tx * pauli.s0sz +1j* p.alphax * pauli.sysz
def hopy(site1, site2, p):
return -p.ty * pauli.s0sz -1j * p.alphay * pauli.sxsz
lat = kwant.lattice.square(a,norbs=4)
syst = kwant.Builder()
syst[(lat(x, y) for x in range(barrierpos)
for y in range(W))]=onsite_barrier
syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy
sym_left = kwant.TranslationalSymmetry((-a, 0))
lead0 = kwant.Builder(sym_left, conservation_law=-pauli.s0sz)
lead0[(lat(0, j) for j in range(W))]=onsite_normal
lead0[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
lead0[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy
sym_right = kwant.TranslationalSymmetry((a, 0))
lead1 = kwant.Builder(sym_right)
lead1[(lat(0, j) for j in range(W))]=onsite_sc
lead1[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
lead1[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy
syst.attach_lead(lead0)
syst.attach_lead(lead1)
syst = syst.finalized()
return syst

Hi,
this is a follow up question to a previous (hijacked?) thread. For completeness I will provide all details, so that no cross referencing is necessary.
I want to compute the electrical conductance of a N-S interface embedded on a hexagonal lattice. The final result should look something like the picture I have attached. The BdG equation is specified in the .pdf.
At the moment I have both particle-hole and spin-up/down degrees of freedom, because later I want to add magnetic fields and such. To attack the problem I have tried to follow the steps outlined in tutorial 2.6 with the addition that I use 4 orbitals and a hexagonal lattice. The conductance should still be given as N - R_{ee,\uparrow} - R_{ee,\downarrow} + R_{he,\uparrow} + R_{he,\downarrow}. Here e, h, \uparrow, and \downarrow refers to electron, hole, spin-up, and spin-down respectively.
In the documentation it says that smatrix.transmission((i, a), (j, b)) gives the transmission from block b of lead j to block a of lead i. I would therefore expect smatrix.transmission((0, 2), (0, 0)) to give me the transmission of an incident electron of spin up in lead 0 to a reflected hole of spin up in lead 0. With this in mind I thought I could write
data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0]-
smatrix.transmission((0, 0), (0, 0)) -
smatrix.transmission((0, 1), (0, 0)) +
smatrix.transmission((0, 2), (0, 0)) +
smatrix.transmission((0,3), (0, 0)))
to compute the conductance. However, doing this gives me an out of bounds error.
As a last resort I tried to write
data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] -
smatrix.transmission((0, 0), (0, 0)) +
smatrix.transmission((0, 1), (0, 0)))
as given in the tutorial 2.6. However, this gives me that the conductance is zero.
My question is then, how do I correctly interpret the smatrix.transmission function in my example so that I can compute the correct conductance?
I have attached my code in the .txt file.
Best,
Martin

Dear Colleagues, (sorry, this is a repost of the same question to a
new thread)
I have a taks to compute local density of states at a single point
for a large homogeneous sample with narrow STS tip potential and in
magnetic field. I have tried several methods in Kwant to do this:
1) Naive diagonalization of H for system with no leads and direct
calculation of LDOS. It works, but not for large systems, so, I see
strong finite-size effects.
2) Sparse Lanczos method with several eigenvalues/functions extracted
near each point in the grid of energy values.
3) KPM method. I am not happy with the energy resolution I can get
with it.
4) Take a smaller main sample (with the size of the tip potential),
but avoid the finite-size effects by attaching a lot of leads pointing
in , e.g. 6 directions. Then, use the kwant ldos function to scan
LDOS over a grid of energies. Use magnetic_gauge to implement B field
in both sample and scattering region.
I see the last possibility as the quickest, but it gives me
something that seems wrong to me.
Does kwant ldos function work correctly if all the leads are insulating
for a given energy? Does it correctly catch the localized states in the
scattering region?
Thank you for any feedback,
Sergey
Below is a copy of a reply by Pablo Piskunow
Dear Sergey,
In terms of scaling, the most convenient is the KPM method, since it scales linearly with the
inverse energy resolution (that is, the number of moments), and linearly with the system size.
Could you clarify what is the issue with the approach 3)?
I am not sure what is the energy resolution that you want to achieve. You can conveniently set
this value when creating a `kwant.kpm.SpectralDensity` instance via the arguments
`energy_resolution` or (exclusive) `num_moments`. The latter is related to the energy resoution
by `\delta \pi / N`, where `\delta` is the full width of the spectrum, and `N` the size of your system.
In the case of a 2D, that will be similar to `L^2`, where `L` is the linear size.
Furthermore, you could progressively increase the energy resolution by adding moments:
```
spectrum = kwant.kpm.SpectralDensity(...)
spectrum.add_moments(...)
```
If you want to resolve individual energy levels then the approach 1) and 3) will have the same scaling,
however, the KPM method will not give you eigenvalues.
Finally, since you want to get the LDOS at a single site or a small set of sites near the STS tip, you should
use the `kwant.kpm.LocalVectors` vector factory, setting the `where` argument to the desired spot.
```
def center(site):
return site.tag == (0, 0)
vector_factory = kwant.kpm.LocalVectors(syst, where=center)
spectrum = kwant.kpm.SpectralDensity(syst, num_moments=100, num_vectors=None, vector_factory=vector_factory)
```
Note that `num_vectors=None` will let the `kwant.kpm` module to figure out how many vectors does the vector factory produce.
Otherwise, it should be at most the*total number of orbitals* defined by the `where` function, that is, sites times orbitals per site.
As you can see, the `kwant.kpm` module is the most suited for this problem, specially when the sample is very large.
I hope this helps.
Best regards,
Pablo
n the meantime, you can create a `spectrum` instance, and evaluate the
KPM expansion at any energy,
for example by
```
spectrum = kwant.kpm.SpectralDensity(...)
energy_array = np.linspace(0, 1, 200)
density_array = spectrum(energy_array)
```
Another note, if the Hamiltonian is huge, you will save some overhead by
passing the bounds of the spectrum explicitly:
`bounds=(e_min, e_max)`. These values must be strictly below and above
the edges of the spectrum, otherwise
the KPM expansion diverges.
ah, I forgot.
The automatic energy points (`spectrum.energies`) are fixed once you fix
the bounds of the spectrum.
So if for a range of magnetic fields you pass the same bounds, then the
energies will be at the same
points and so the densities will be easier to plot or compare.
The width of the spectrum might change for different values of the
magnetic field; and even if not, there is
small randomness associated with finding the bounds, that you can get
rid of by passing a seed to the
random number generator like `rng=0`. If you choose to provide the
bounds explicitly, you can disregard
this random number generator stuff.

Hello,
I have recently started using KWANT. I hope this is the correct way to ask questions.
I'm interested in reproducing Fig. 4 in PRL 97, 067007 (2006) - Specular Andreev Reflection in Graphene. This has been done in KWANT in the PhD thesis of Tibor Sekera - "Quantum transport of fermions in honeycomb lattices and cold atomic systems", chapter 4. The problem is similar to tutorial 2.6.
In my problem, I have a tight-binding Hamiltonian living on a hexagonal lattice, from which I can construct a discrete BdG equation (See attachment). In addition, to electron-hole degrees of freedom I'm also interested in including spin because later I will add magnetic fields and such.
I consider a normal metal lead attached to a superconducing scattering region with zero barrier at the interface. Since I have both electron-hole and spin degrees of freedom my BdG equation is 4x4. To account for electron-hole degrees of freedom I have specified that at each site there should be two orbitals. In the tutorial 2.6 it is emphasized that it is important to specify the correct conservation law in the N lead when we want to calculate the conductance. My question is how should the conservation law (and particle hole symmetry) in the normal-metal lead look like for my case?
I have tried the standard expression L_0 = kwant.Builder(sym,conservation_law=s_z) where s_z is the Pauli matrix [[1,0],[0,-1]]. In this case I obtain a dimension mismatch error.
I have also tried a generalized 4x4 matrix version of this
L_0 = kwant.Builder(sym,conservation_law=tinyarray.array([[1,0,0,0],[0,-1,0,0],[0,0,1,0],[0,0,0,-1]]))
Or
L_0 = kwant.Builder(sym,conservation_law=tinyarray.array([[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]])).
In these cases I obtain the error "Single `onsite` matrix of shape (4, 4) provided but there are 2 orbitals per site in the system".
I think my issue is that I do not fully understand what the Conservation law actually does. Initially I thought the conservation law was a matrix C such that C.H.C^-1 should return a block diagonal matrix, where the eigenvalues of C refers to each block. Is this not an accurate interpretation?
I have attached a PDF with the BdG equation, and my python code as a .py and .txt file. The content of the .py and .txt file are identical, which format do you prefer in the future?
I hope this is OK, and not too overwhelming.
I also tried to subscribe (with this e-mail) to the mailing list, but I'm not sure that it worked.
Yours sincerely,
Martin F Jakobsen

I know there is Tbmodels package can deal with wannier90 output tight binding files. But I think TBmodels group all the orbitals into one site.
I was trying to do things that can assign hoppings and on-site energies in "Real" systems that have differnet sites but with different orbitals per site.
The errors are the following when I try to use kwant_model.hamiltonian_submatrix module
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/builder.py in hamiltonian(self, i, j, params, *args)
1885 try:
-> 1886 value = value(site, *args)
1887 except Exception as exc:
/Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/wraparound.py in f(*in_args)
154 if callable(val):
--> 155 acc = acc + val(*out_args)
156 else:
ValueError: operands could not be broadcast together with shapes (6,6) (10,10)
The above exception was the direct cause of the following exception:
UserCodeError Traceback (most recent call last)
<ipython-input-34-1b24fb1489bf> in <module>()
3 params={key: val for key, val in zip(['k_x', 'k_y', 'k_z'], 2 * np.pi * np.array(k_list[0]))}
4 )
----> 5 ) for k in k_list]
<ipython-input-34-1b24fb1489bf> in <listcomp>(.0)
3 params={key: val for key, val in zip(['k_x', 'k_y', 'k_z'], 2 * np.pi * np.array(k_list[0]))}
4 )
----> 5 ) for k in k_list]
/Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/_common.py in inner(*args, **kwargs)
70 if sig.bind(*args, **kwargs).arguments.get(parameter_name):
71 warn()
---> 72 return f(*args, **kwargs)
73
74 return inner
kwant/_system.pyx in kwant._system.hamiltonian_submatrix()
/Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/builder.py in hamiltonian(self, i, j, params, *args)
1892 ', '.join(map('"{}"'.format, missing)))
1893 raise TypeError(''.join(msg))
-> 1894 _raise_user_error(exc, value)
1895 else:
1896 edge_id = self.graph.first_edge_id(i, j)
/Users/x0k/opt/anaconda3/envs/kwant/lib/python3.6/site-packages/kwant/builder.py in _raise_user_error(exc, func)
1814 msg = ('Error occurred in user-supplied value function "{0}".\n'
1815 'See the upper part of the above backtrace for more information.')
-> 1816 raise UserCodeError(msg.format(func.__name__)) from exc
1817
1818
UserCodeError: Error occurred in user-supplied value function "f".
See the upper part of the above backtrace for more information.

Dear Kwant Developers,
I would like to know the unit of the current density shown in the current plot. Does it always equal "unit of current/relwidth" for a 2D system? I understand that abswidth can replace relwidth. If I do not specify any abswidth, will the relwidth be replaced? Finally, the relwidth is defined as , as a fraction of the length of the longest side of the bounding box. What is the definition of the bounding box? Is it the area in which the current is plotted? Thank you in advance for your help.
Regards,
KS Chan
Disclaimer: This email (including any attachments) is for the use of the intended recipient only and may contain confidential information and/or copyright material. If you are not the intended recipient, please notify the sender immediately and delete this email and all copies from your system. Any unauthorized use, disclosure, reproduction, copying, distribution, or other form of unauthorized dissemination of the contents is expressly prohibited.

Dear all,
I am trying a graphene bulk with magnetic field, and I use
spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None,
vector_factory=s_factory,
rng=0)
energies, densities = spectrum()
to calculate the dos and compare with some previous studies. My definition
of the hopping is:
def Hop_magnetic(site1,site2):
B_magneic=200. # units Tesla
a0=0.142 # nm
x1,y1=site1.pos
x2,y2=site2.pos
xy=0.5*(x1+x2)*sqrt(3)*a0*(y1-y2)*sqrt(3)*a0 # the units is nm**2
phb=1j*2.*pi*B_magneic/4.135667 # phi0=h/e=4.135667*e-15 V*s
return t*exp(xy*phb*0.001)
I found that the effect of the magnetic field is too small and I can not
obtain the correct results.
Could you give me some suggestions?
Best regards
Khani
My code is pasted:
import kwant
from matplotlib import pyplot
from numpy import sqrt,pi,exp
import numpy as np
def make_system(r=8, t=1, tp=-0.1):
lat = kwant.lattice.honeycomb(norbs=1)
a, b = lat.sublattices
def circle(pos):
x, y = pos
return x**2 + y**2 < 100**2 #-100<x<100 and -100<y<100
def Hop_magnetic(site1,site2):
B_magneic=200. # the unit of B_magneic is Tesla=V*s/m**2
a0=0.142 # nm
x1,y1=site1.pos
x2,y2=site2.pos
xy=0.5*(x1+x2)*sqrt(3)*a0*(y1-y2)*sqrt(3)*a0 # the units change
it to nm**2
phb=1j*2.*pi*B_magneic/4.135667 # phi0=h/e=4.135667*e-15 V*s,
the unit of B_magneic is Tesla=V*s/m**2
return t*exp(xy*phb*0.001)
def Nearest_Hop(site1,site2):
return Hop_magnetic(site1,site2)
syst = kwant.Builder()
syst[lat.shape(circle, (0, 0))] = 0
syst[lat.neighbors()] = Nearest_Hop
syst.eradicate_dangling()
return lat, syst.finalized()
lat, fsyst = make_system()
#kwant.plot(fsyst)
where = lambda s: np.linalg.norm(s.pos) < 1
s_factory = kwant.kpm.LocalVectors(fsyst, where)
spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None,
vector_factory=s_factory,
rng=0)
energies, densities = spectrum()
pyplot.figure()
pyplot.plot(energies, densities)
pyplot.xlabel("energy [t]")
pyplot.ylabel("DOS")
pyplot.show()