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

Hello everyone
I want to use different lattices structure (honeycomb and diamond) for
the scattering region. if I have a honeycomb structure, is it possible
in Kwant to add sites from the diamond structure to the system and add
hoppings from these sites to the sites from the honeycomb lattice?
I have already read Section 2.11
(https://kwant-project.org/doc/1/tutorial/faq#how-to-use-different-lattices-…)
in the documentation but this is challenging for me and I have no idea
how to solve it.
Thanks

Dear Kwant developers and users,
I write concerning a doubt that I encountered
trying to study by Kwant a model for the quantum Hall/quantum spin Hall effect.
In particular, I am focusing on the two-dimensionakl BHZ model, first
proposed in the famous paper
https://science.sciencemag.org/content/314/5806/1757
Notice that first I keep only one of the two time-reversal-connected blocks
forming the quantum spin Hall system (when the mass parameter M is positive);
it is known that each block describes a quantum Hall system (again for
M >0, otherwise a trivial insulator occurs).
I attach below my defining code (included in the % lines)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
import kwant
from kwant.digest import uniform # a (deterministic) pseudorandom
number generator
import math
from matplotlib import pyplot as plt
import cmath
import scipy.sparse.linalg as sla
import scipy.linalg as la
import numpy
from numpy import linalg as LA
# For plotting
from matplotlib import pyplot
# For matrix support
import tinyarray
#import progressbar
#-----------------------------------------------------
# We define the variables and the constants
a = 1 #lattice step
#t = 1 #along x and y
#-----------------------------------------------------
nullmatrix = tinyarray.array([[0, 0], [0, 0]])
matid = tinyarray.array([[1, 0], [0, 1]])
sigma_0 = tinyarray.array([[1, 0], [0, -1]])
sigma_xp = tinyarray.array([[1, -1j], [-1j, -1]])
sigma_yp = tinyarray.array([[1, 1/2], [-1/2, -1]])
-----------------------------------------------------------------------------------------------------------------
def sistema_BHZ_2e(Lx, Ly, M) :
syst = kwant.Builder()
lat = kwant.lattice.square(norbs=2)
for i in range(0,Lx):
#syst[lat(i, 0), lat(i, L - 1)] = sigma_yp
for j in range(0,Ly):
# On-site Hamiltonian
syst[lat(i, j)] = (M - 4)*sigma_0
# Hopping in y-direction
if i > 0:
syst[lat(i, j), lat(i - 1, j)] = sigma_xp
# Hopping in z-direction
if j > 0:
syst[lat(i, j), lat(i, j - 1)] = sigma_yp
return syst
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Correspondingly to what described above, I would expect to see, when M
>0, edge modes.
Indeed, I obtain them calculating the local densities via the KPM
(approximate) method working in Kwant:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Plot several local density of states maps in different subplots
def plot_ldos(syst, densities, file_name=None):
step = 0.1
fig, axes = plt.subplots(1, len(densities), figsize=(7*len(densities), 7))
for ax, (title, rho) in zip(axes, densities):
kwant.plotter.density(syst,rho.real, colorbar=True, ax=ax,
relwidth = step)
ax.set_title(title)
ax.set(adjustable='box', aspect='equal')
plt.show()
plt.clf()
--------------------------------------------------------------------------------------------------------------
M = 2
L = 16
Lx = L
Ly = L
syst = sistema_BHZ_2e(Lx, Ly, M)
syst = syst.finalized()
hamat = syst.hamiltonian_submatrix(sparse = False)
evals, evecs = la.eigh(hamat)
kwant_op = kwant.operator.Density(syst, sum=False)
local_dos = kwant.kpm.SpectralDensity(syst, operator=kwant_op)
print("first negative energy =", evals[int(Lx*Ly-1)])
print("first_positive_energy =", evals[int(Lx*Ly)])
first_negative_energy_ldos = local_dos(energy=evals[int(Lx*Ly-1)])
first_positive_energy_ldos = local_dos(energy=evals[int(Lx*Ly)])
plot_ldos(syst, [('first negative energy',
first_negative_energy_ldos),('first positive energy',
first_positive_energy_ldos)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I attach below a typical output (KPM.pdf).
However, when I calculate the same densities not using the KPM approach
(then via the exact function kwant.operator.Density(syst) only ) ,
that means by the code:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = 2
L = 16
Lx = L
Ly = L
syst = sistema_BHZ_2e(Lx, Ly, M)
syst = syst.finalized()
hamat = syst.hamiltonian_submatrix(sparse = False)
evals, evecs = la.eigh(hamat)
kwant_op = kwant.operator.Density(syst, sum=False)
print("first negative energy =", evals[int(Lx*Ly-1)])
print("first_positive_energy =", evals[int(Lx*Ly)])
first_negative_energy_ldos = kwant_op(evecs[int(Lx*Ly-1)])
first_positive_energy_ldos = kwant_op(evecs[int(Lx*Ly)])
plot_ldos(syst, [('first negative energy',
first_negative_energy_ldos),('first positive energy',
first_positive_energy_ldos)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I obtain a totally different density profile, not displaying
edge localization (see No_KPM.pdf attached below).
I do not understand the reason of such a mismatch.
Can you help me ?
Thank you very much and best,
L.

Dear all,
I want to discrete the continuous Hamiltonians in PRL 111, 156402 (2013)
and calculate the band structure of the ribbon. So I change the parameters
of the code in "Tutorial 2.10. Discretizing continuous Hamiltonians". But I
can not obtain the same band structure as shown in Fig.3(b) of the PRL
paper. Could you please help me to check my attached code? Thanks!
Yours sincerely
Hosein Khani
###########################
import kwant
import kwant.continuum
import scipy.sparse.linalg
import scipy.linalg
import numpy as np
# For plotting
import matplotlib as mpl
from matplotlib import pyplot as plt
def qsh_system(a=20, L=4000, W=5000):
hamiltonian = """
+ C * identity(4) + M * kron(sigma_0, sigma_z)
- B * (k_x**2 + k_y**2) * kron(sigma_0, sigma_z)
- D * (k_x**2 + k_y**2) * kron(sigma_0, sigma_0)
+ A * k_x * kron(sigma_z, sigma_x)
- A * k_y * kron(sigma_0, sigma_y)
"""
template = kwant.continuum.discretize(hamiltonian, grid=a)
def shape(site):
(x, y) = site.pos
return (0 <= y < W and 0 <= x < L)
def lead_shape(site):
(x, y) = site.pos
return (0 <= y < W)
syst = kwant.Builder()
syst.fill(template, shape, (0, 0))
lead = kwant.Builder(kwant.TranslationalSymmetry([-a, 0]))
lead.fill(template, lead_shape, (0, 0))
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
syst = syst.finalized()
return syst
def analyze_qsh():
#######the follwing parameters are given to make sure the hamiltonian is
the same as
####### PRL 111, 156402 (2013)
AP=0.02851
BP=(0.4381-0.2081)/2.
DP=(0.4381+0.2081)/2.
MP=(-0.19808+0.19153)/2.
CP=(-0.19808-0.19153)/2.
params = dict(A=AP, B=BP, D=DP, M=MP, C=CP)
syst = qsh_system()
kwant.plotter.bands(syst.leads[0], params=params,
momenta=np.linspace(-0.3, 0.3, 201), show=False)
plt.grid()
plt.xlim(-.3, 0.3)
plt.ylim(-0.2, -0.19)
plt.xlabel('momentum [1/A]')
plt.ylabel('energy [eV]')
plt.show()
# get scattering wave functions at E=0
wf = kwant.wave_function(syst, energy=0, params=params)
# prepare density operators
sigma_z = np.array([[1, 0], [0, -1]])
prob_density = kwant.operator.Density(syst, np.eye(4))
spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2)))
# calculate expectation values and plot them
wf_sqr = sum(prob_density(psi) for psi in wf(0)) # states from left
lead
rho_sz = sum(spin_density(psi) for psi in wf(0)) # states from left
lead
kwant.plotter.map(syst, wf_sqr, show=False)
plt.show()
def main():
analyze_qsh()
# Call the main function if the script gets executed (as opposed to
imported).
# See <http://docs.python.org/library/__main__.html>.
if __name__ == '__main__':
main()

Hi everyone!
Together with the other Kwant authors I am planning a grant
application. One activity for which we consider requesting support is
a series of internships—visits to one of our groups with the goal of
developing either a Kwant feature or a related project. These
internships would last for perhaps a couple of months, with the idea
that by the end of the visit, they result in something either
implemented or being close to completion.
We think that this work organization would benefit everyone involved,
and we'm very much interested to explore this idea.
Still, before we commit to that, we would be very much interested to
hear everyone's opinion: Would you or someone you know want to go for
such an internship? What do you think of this idea overall? Any idea
of what project could be useful?
Cheers,
Anton

Dear Kwant developers,
What are the units of the output of the current density operator? The units of the density operator is per area per energy, so the units of current are (e/h)* per area * energy?
Thanks a lot.
Bests,
Marc

Dear all,
Now I’m studying some transport property based on multilayered system, with nonmagnetic layer/ferromagnetic layer (NM/FM) stack.
May I ask when creating a tight binding system, and setting hoppings by the following function:
syst[(lat(1, 0), lat(1, 1))] = t,
how to associate the hopping with a realistic system? Such as hoppings between FM and FM atoms, or between FM and NM atoms? I found in many literatures people just simply put t=1, but I assume the value should depend on the type of atoms.
My another question is, is it possible to define more than one type of hopping between two sites?
For the tight-binding approximation of p electrons, there are two types of hopping integrals: Vppσ, and Vppπ, so I’m wondering if it’s possible to simulate the effect of both hoppings?
Thanks a lot for reading my questions!
Cheers,
Yu Li

Hello Kwant community,
I write since I am trying to simulate a Kitaev chain with disorder,
here distributed uniformly in the range [1,-1],
and I am experiencing troubles to implement the disorder itself.
In the absence of disorder,
the function that defines the tight-binding (open) Hamiltonian can be written
e.g. as follows:
def kitaev(L ,t, mu, Delta) :
syst = kwant.Builder()
lat = kwant.lattice.chain(norbs=2)
syst[lat(0)] = - (mu/2)*pauli3
for i in range(1,L):
syst[lat(i)] = -(mu/2)*pauli3
syst[lat(i), lat(i - 1)] = -(t/2)*pauli3 + 1j*(Delta/2)*pauli2
return syst
where the Pauli matrices have been previously defined, t is the hopping,
mu is the chemical potential, and Delta is the gap.
Now, I want to add disorder, proportional to W and uniformly
distributed in [-1,1]. The tutorial online for Kwant states that the
function kwant.digest.uniform(input, salt='') creates a random offset
in the interval [0,1]. Therefore I thought that the Hamiltonian could
be changed e.g. as follows:
def kitaev(L ,t, mu, Delta, W) :
syst = kwant.Builder()
lat = kwant.lattice.chain(norbs=2)
syst[lat(0)] = -(1/2)*(mu + W*(2*onsite - 1))*pauli3
for i in range(1,L):
syst[lat(i)] = -(1/2)*(mu + W*(2*onsite - 1))*pauli3
syst[lat(i), lat(i - 1)] = -(t/2)*pauli3 + 1j*(Delta/2)*pauli2
return syst
after that the onsite function has been defined as
def onsite(site, phi, salt):
return uniform(repr(site), salt)
(I am still trying to follow the tutorial).
However, when trying to calculate eigenvalues and eigenvectors,
via the code
L = 50
W = 0
mub = 3
systb = kitaev(L ,1, mub, 1, W)
systb = systb.finalized()
hamatb = systb.hamiltonian_submatrix(sparse = False)
evalsb, evecsb = la.eigh(hamatb)
I encountered the problem
-----------------------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-10-12a574eefcd0> in <module>
3
4 mub = 3
----> 5 systb = kitaev(L ,1, mub, 1, W)
6 systb = systb.finalized()
7 hamatb = systb.hamiltonian_submatrix(sparse = False)
<ipython-input-9-661d5ddb539a> in kitaev(L, t, mu, Delta, W)
4 lat = kwant.lattice.chain(norbs=2)
5
----> 6 syst[lat(0)] = -(1/2)*(mu + W*(2*onsite - 1))*pauli3
7
8 for i in range(1,L):
TypeError: unsupported operand type(s) for *: 'int' and 'function'
—————————————————————————————------
that instead does no occur in the absence of disorder (then when the
first definition above for the Kitaev chain is assumed).
Probably I did not understand properly how to use the
function uniform(repr(site), salt).
Can you help me ?
Thank you very much in advance and best regards
L. L.