Dear Johnny,

The error message you're seeing isn't lying. You probably have a lead
with zigzag boundary conditions that has a bad singularity at E=0.
Choosing even E=1e-4 would resolve the problem. And yes, you're right
about the spin density (but of course if you don't add any
spin-dependent terms, there won't be any spin physics happening.

Best regards,
Anton Akhmerov

On Sun, Nov 29, 2015 at 11:28 AM, T.C. Wu <tcwu@connect.ust.hk> wrote:
> 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
> #execfile('C:/Users/wtc/Google Drive/Linnian/plot_graphene.py')
> 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
>
>     #### Define the leads. ####
>     # left lead
>     sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
>
>     def lead0_shape(pos):
>         x, y = pos
>         return (-0.4 * r < y < 0.4 * r)
>
>     lead0 = kwant.Builder(sym0)
>     lead0[graphene.shape(lead0_shape, (0, 0))] = 0
>     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)))
>
>     def lead1_shape(pos):
>         v = pos[1] * sin_30 - pos[0] * cos_30
>         return (-0.4 * r < v < 0.4 * r)
>
>     lead1 = kwant.Builder(sym1)
>     lead1[graphene.shape(lead1_shape, (0, 0))] = 0
>     lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] =
> hopping2
>
>     return sys, [lead0, lead1]
>
>
> 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()
>
>
> def plot_bandstructure(flead, momenta):
>     bands = kwant.physics.Bands(flead)
>     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()
> def density(sys, energy, args, lead_nr):
>     wf = kwant.wave_function(sys, energy, args)
>     return (abs(wf(lead_nr))**2).sum(axis=0)
>
> def main():
>     pot = 0.1
>     energy = 0.0
>     sys, leads = make_system(pot=pot)
>
>     # 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.
>     for lead in leads:
>         sys.attach_lead(lead)
>
>     # Then, plot the system with leads.
>     kwant.plot(sys, site_color=family_colors, site_lw=0.1,
>                lead_site_lw=0, colorbar=False)
>
>     # 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_bandstructure(sys.leads[0], momenta)
>
>
>     # 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, ""])
>         conductances.append(smatrix.transmission(0, 1)) # from lead 1 to
> lead 0
>     pyplot.plot(reciprocal_phis, conductances)
>     pyplot.show()
>     # Plot conductance.
>     # energies = [-2 * pot + 4. / 50. * pot * i for i in xrange(51)]
>     # plot_conductance(sys, 0)
>
>
> # 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()
>