Dear all,

I am trying to use kwant to reproduce the QHE in graphene. I combine the example in qhe.py and graphene.py provided in the official website(see the code below).However, I got the following error when I plot out the wavefunction at energy = 0:
RuntimeError: Numbers of left- and right-propagating modes differ, possibly due to a numerical instability.

Moreover, if I don't choose energy = 0, say energy = 0.1, then the error disappear. However, I cannot recover the plateaus in the plot of conductance.

Lastly, it is correct to change t -> t*sigma_0 in order to extract the distribution of spin(i.e. choosing the odd /even elements in the wavefunction by [1::2]/[0::2]) in this case?

Thank you.

Best,
Johnny

Code:
# Tutorial 2.5. Beyond square lattices: graphene
# ==============================================
#
# Physics background
# ------------------
#  Transport through a graphene quantum dot with a pn-junction
#
# Kwant features highlighted
# --------------------------
#  - Application of all the aspects of tutorials 1-3 to a more complicated
#    lattice, namely graphene
from __future__ import division  # so that 1/2 == 0.5, and not 0
from math import pi, sqrt, tanh
import random
import kwant
import numpy as np
# For computing eigenvalues
import scipy.sparse.linalg as sla

# For plotting
from matplotlib import pyplot

# Define the graphene lattice
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
[(0, 0), (0, 1 / sqrt(3))])
a, b = graphene.sublattices
Bs = np.linspace(4,50,200)
t = -1

def make_system(r=10, w=2.0, pot=0.1):

#### Define the scattering region. ####
# circular scattering region
def circle(pos):
x, y = pos
return x ** 2 + y ** 2 < r ** 2
#return y > (abs(x)*np.tan(np.pi/3) -size) and y < size
sys = kwant.Builder()

# w: width and pot: potential maximum of the p-n junction
def potential(site):
(x, y) = site.pos
d = y * cos_30 + x * sin_30
return pot * tanh(d / w)

def hopping2(site1, site2, phi = 0,salt=0):
xt, yt = site1.pos
xs, ys = site2.pos
return -t * np.exp(-1j * phi * (xt-xs) *(yt+ys))

sys[graphene.shape(circle, (0, 0))] = 0

# specify the hoppings of the graphene lattice in the
# format expected by builder.HoppingKind
hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = hopping2

# Modify the scattering region
# del sys[a(0, 0)]
# del sys[a(5, 5)]
# del sys[a(5,5)]
# sys[a(-2, 1), b(2, 2)] = -1

sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))

x, y = pos
return (-0.4 * r < y < 0.4 * r)

lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = hopping2

# The second lead, going to the top right
sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))

v = pos[1] * sin_30 - pos[0] * cos_30
return (-0.4 * r < v < 0.4 * r)

lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = hopping2

def compute_evs(sys):
# Compute some eigenvalues of the closed system
sparse_mat = sys.hamiltonian_submatrix(sparse=True)

evs = sla.eigs(sparse_mat, 2)[0]
print evs.real

def plot_conductance(sys, energy):
# Compute transmission as a function of energy
data = []
for phi in Bs:
smatrix = kwant.smatrix(sys, energy,args =[phi,""])
data.append(smatrix.transmission(0, 1))

pyplot.figure()
pyplot.plot(Bs, data)
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()

energies = [bands(k) for k in momenta]

pyplot.figure()
pyplot.plot(momenta, energies)
pyplot.xlabel("momentum [(lattice constant)^-1]")
pyplot.ylabel("energy [t]")
pyplot.show()
wf = kwant.wave_function(sys, energy, args)

def main():
pot = 0.1
energy = 0.0

# To highlight the two sublattices of graphene, we plot one with
# a filled, and the other one with an open circle:
def family_colors(site):
return 0 if site.family == a else 1

# Plot the closed system without leads.
kwant.plot(sys, site_color=family_colors, site_lw=0.1, colorbar=False)

# Compute some eigenvalues.
compute_evs(sys.finalized())

# Attach the leads to the system.

# Then, plot the system with leads.
kwant.plot(sys, site_color=family_colors, site_lw=0.1,

# Finalize the system.
sys = sys.finalized()

# Compute the band structure of lead 0.
momenta = [-4*pi + 0.02 * 4*pi * i for i in xrange(101)]

# plot wf
d = density(sys, energy, [1/40.0, ""], 0)
kwant.plotter.map(sys, d)
reciprocal_phis = np.linspace(4, 50, 200)
conductances = []
for phi in 1 / reciprocal_phis:
smatrix = kwant.smatrix(sys, energy, args=[phi, ""])