Hello Kwant community,
I am using the KPM method to calculate conductivities in the simplest QHE system -- square lattice nearest hopping gas. The goal is to see well-quantized plateaus up to n=5.
I use local vectors at the center of the system. But it seems not available at all as far as I've tried increasing system size (like 80*80) and number of moments (like 500). I doubt if it really needs so demanding computational resources in order to get n=5. Hopefully I missed something basic here.
Thank you!
from cmath import exp
import numpy as np
import time
import sys
import kwant
#from matplotlib import pyplot
from matplotlib import pyplot as plt
L_sys = [50, 50]
E_F=0.4
def make_system(a=1, t=1):
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
def fluxx(site1, site2, Bvec):
pos = [(site1.pos[i]+site2.pos[i]-L_sys[i])/2.0 for i in range(2)]
return exp(-1j * Bvec * pos[1] )
def fluxy(site1, site2, Bvec):
pos = [(site1.pos[i]+site2.pos[i]-L_sys[i])/2.0 for i in range(2)]
return 1 #exp(-1j * (-Bvec * pos[0]))
def hopx(site1, site2, Bvec):
# The magnetic field is controlled by the parameter Bvec
return fluxx(site1, site2, Bvec) * t
def hopy(site1, site2, Bvec):
# The magnetic field is controlled by the parameter Bvec
return fluxy(site1, site2, Bvec) * t
def onsite(site):
return 4*t
#### Define the scattering region. ####
syst[(lat(x, y) for x in range(L_sys[0]) for y in range(L_sys[1]))] = onsite
# hoppings in x-direction
syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
# hoppings in y-directions
syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy
# Finalize the system.
syst = syst.finalized()
return lat, syst
def cond(lat, syst, random_vecs_flag=False):
center_pos = [x/2 for x in L_sys]
center0 = lambda s: np.linalg.norm(s.pos-center_pos) < 2
num_mmts = 800;
if random_vecs_flag:
center_vecs = kwant.kpm.RandomVectors(syst, where=center0)
one_vec = next(center_vecs)
num_vecs = 10 # len(center_vecs)
else:
center_vecs = np.array(list(kwant.kpm.LocalVectors(syst, where=center0)))
one_vec = center_vecs[0]
num_vecs = len(center_vecs)
area_per_orb = np.abs(np.cross(*lat.prim_vecs))
norm = area_per_orb * np.linalg.norm(one_vec) ** 2
Binvfields = np.arange(3.01, 26.2, 1)
params = dict(Bvec=0)
dataxx = []
dataxy = []
for Binv in Binvfields:
params['Bvec'] = 1/Binv
cond_xx = kwant.kpm.conductivity(syst, params=params, alpha='x', beta='x', mean=True,
num_moments=num_mmts, num_vectors=num_vecs, vector_factory=center_vecs)
cond_xy = kwant.kpm.conductivity(syst, params=params, alpha='x', beta='y', mean=True,
num_moments=num_mmts, num_vectors=num_vecs, vector_factory=center_vecs)
dataxx.append(5*cond_xx(mu=E_F, temperature=0.0)/norm)
dataxy.append(-cond_xy(mu=E_F, temperature=0.0)/norm)
print(dataxx)
print(dataxy)
plot_curves(
[
(r'$\sigma_{xx}\times5$',(Binvfields,dataxx)),
(r'$\sigma_{xy}$',(Binvfields,dataxy)),
]
)
def plot_curves(labels_to_data):
plt.figure(figsize=(12,10))
for label, (x, y) in labels_to_data:
plt.plot(x, y, label=label, linewidth=2)
plt.legend(loc=2, framealpha=0.5)
plt.xlabel("1/B")
plt.ylabel(r'$\sigma [e^2/h]$')
plt.show()
plt.clf()
def main():
lat, syst = make_system()
cond(lat, syst)
if __name__ == '__main__':
main()

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()

Dear author：
I am a PHD student and my research is using pizeoelectric potential
regulate spin and valley transport. I am quite intristed about Kwant and I
have published an article on Nona Energy using Kwant
I wonder could you please help me to solve the folowing two
questions, the first one is how to calculate a one-dimensional quantum
transport system, can you give a case? The second one is how to use Kwant
to caculate heat transport system like electric current, heat current,
temperature-dependent electrical conductance , the thermal conductivity,
the Seebeck thermopower, and the thermoelectric efficiency. Could you give
me an example? Thank you.
Best wishes
Ruhao Liu

Hello,
I am trying to extract Green's function between two arbitrary sites in a system. I have noticed from one of the previous threads that one way to do it is to attach virtual leads with zero self-energy on those specific sites. However, I am having trouble with the suggested function. It is telling me that whichever site I want to attach is not in the scattering region. The following is the code:
import kwant
from types import SimpleNamespace
from copy import copy
import tinyarray
import numpy as np
import csv
import ipywidgets
#import holoviews as hv
#hv.extension('bokeh')
from scipy import sparse
from matplotlib import pyplot
from scipy.optimize import minimize
sigma_x = tinyarray.array([[0, 1], [1, 0]])
sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
sigma_z = tinyarray.array([[1, 0], [0, -1]])
sigma_0 = tinyarray.array([[1, 0], [0, 1]])
tau_x = tinyarray.array([[0, 1], [1, 0]])
tau_y = tinyarray.array([[0, -1j], [1j, 0]])
tau_z = tinyarray.array([[1, 0], [0, -1]])
tau_0 = tinyarray.array([[1, 0], [0, 1]])
# mu=0.5,
# t=-1.,delta=0.5,V_N=1.25, phi=0, L=3
def make_system_ex3(L=10):
lat = kwant.lattice.chain(norbs=2)
syst = kwant.Builder()
#### Define the scattering region. ####
def onsite_sc(site,t,mu,delta):
return (2.*t-mu )*tau_z + delta*tau_x
def onsite_N(site,t,mu,V_N):
return (2.*t-mu + V_N)*tau_z
syst[(lat(x) for x in range(1, L))] = onsite_sc
#superconducting onsite hamiltonian
syst[lat(L)] = onsite_N #normal onsite hamiltonian
syst[(lat(x) for x in range(L+1,2*L+1))] = onsite_sc #superconducting onsite hamiltonian
def hop(site1,site2,t,phi):
return -t*np.matmul(np.exp(1j*phi*tau_z),tau_z)
syst[(lat(L-1), lat(L))] = hop
for i in range(1,L-1):
syst[(lat(i), lat(i+1))] = -t*tau_z
for i in range(L,2*L):
syst[(lat(i), lat(i+1))] = -t*tau_z
#### Define the leads. ####
sym_left = kwant.TranslationalSymmetry([-1])
lead0 = kwant.Builder(sym_left)
lead0[(lat(0))] = (2.*t-mu)*tau_z + delta*tau_x
lead0[lat.neighbors()] = -t*tau_z
sym_right = kwant.TranslationalSymmetry([1])
lead1 = kwant.Builder(sym_right)
lead1[(lat(2*L+1))] = (2.*t-mu)*tau_z + delta*tau_x
lead1[lat.neighbors()] = -t*tau_z
# #### Attach the leads and return the system. ####
syst.attach_lead(lead0)
syst.attach_lead(lead1)
return syst,lat
----------------------------------------------------------------------------------------------------------------------------
mu=1
t=1.0
delta=0.5
V_N=1
phi=0
L=10
par = dict(t=t,mu=mu,delta=delta,phi=phi,V_N=V_N)
#
syst,lat = make_system_ex3(L=10)
def mount_vlead(sys, vlead_interface, norb):
dim = norb
zero_array = np.zeros((dim, dim), dtype=float)
def selfenergy_func(energy, args=()):
return zero_array
vlead = kwant.builder.SelfEnergyLead(selfenergy_func, vlead_interface,())
sys.leads.append(vlead)
mount_vlead(syst,[L-1], 2) # number of orbitals =4. Here we mount lead number 2
mount_vlead(syst,[L], 2) # number of orbitals =4. Here we mount lead number 3
syst =syst.finalized()
G12=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[L-1],out_leads=[L],\
check_hermiticity=False,params=par).data
G21=kwant.greens_function(syst, energy=-1.8*1j, in_leads=[L],out_leads=[L-1],\
check_hermiticity=False,params=par).data
----------------------------------------------------------------------------------------------------------------------------------
The error message I am getting is
" ValueError: Lead 2 is attached to a site that does not belong to the scattering region:
9"
It would be great if you can suggest to me where I am going wrong.
Thanks,
Sayandip

Dear developer,
I am looking for the calculation of spin Hall conductance in a four-terminal device. I am searching the mailing archives for a few days, but couldn't get anything useful. But, I found a thread where Mr Joseph Weston (https://mail.python.org/archives/list/kwant-discuss@python.org/message/MHYO…) had given some links for this calculation and unfortunately those links are not working. It would be very useful if you provide one example of calculating spin Hall conductance of a four-terminal device.
I shall look forward to your kind reply.
Thanking you.
Sincerely yours,
Sayan Mondal,
IIT Guwahati,