Hello, I was successful in introducing disorder in my Haldane (second nearest neighbours) system. I am using kwant1.3, I have read the new two threads [1] and [2], and I have read [3]. However, I am still unable to compute and plot the local density current for my system in the scattering region. I would really appreciate any help. I attach the code. Once I have it I will attach it, so people can actually use it and solve future problems. Regards, Enrique [1] https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html> [2] https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html> [3] <https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html> https://kwant-project.org/doc/dev/reference/generated/kwant.operator.Current from scipy.sparse import spdiags from math import pi, sqrt, tanh, e import numpy as np import kwant import random import math import tinyarray import scipy.sparse.linalg as sla from matplotlib import pyplot # Define the graphene lattice graphene = kwant.lattice.honeycomb() a, b = graphene.sublattices L=4 W=4 t=1 tt= 0.03 * t phi=1*pi/2 def make_system(pot=0.000): nnn_hoppinga=tt*e**(1j*phi) nnn_hoppingb=tt*e**(1j*phi) #### Define the scattering region. #### def rec(pos): x, y = pos return 0 < x < L and 0 < y < W sys = kwant.Builder() # w: width and pot: potential maximum of the p-n junction def potential(site): (x, y) = site.pos d = y * 2 + x * 3 return 0 sys[graphene.shape(rec, (1, 1))] = potential # specify the hoppings of the graphene lattice nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb #### Define the leads. #### # left lead sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0))) sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead0_shape(pos): x, y = pos return (0 < y < W) lead0 = kwant.Builder(sym0) lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb # right lead sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0))) sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead1_shape(pos): x, y = pos return (0 < y < W) lead1 = kwant.Builder(sym1) lead1[graphene.shape(lead1_shape, (0, 0))] = pot lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb sys.attach_lead(lead1) sys.attach_lead(lead0) for site in sys.leads[0].interface+sys.leads[1].interface: sys[site]=0 return sys def jcurrent(sys): J=kwant.operator.Current(sys) wfs = kwant.wave_function(sys) wf= wfs[0] j=J(wf) kwant.plotter.map(sys, j) def main(): sys = make_system() # Finalize the system. sys = sys.finalized() # Plot local current jcurrent(sys) # Call the main function if the script gets executed (as opposed to imported). # See <http://docs.python.org/library/__main__.html><http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main() El 09/02/2017 a las 2:04, Joseph Weston escribió: Hi again, b) When I asked about computing it directly from Kwant it was because I read in another thread that it was going to be in Kwant1.3, but I was not sure if it was already avalaible. I will "fight" then with the code and will try to use your guidelines Joseph. Definetively doing it with neighbors will be much faster than itering the whole system. Will this new version be ready soon? :-) Well it is already available in the 'development' version of Kwant. If you're using 'conda' to install kwant you can easily install the development version with: conda install -c kwant kwant otherwise you can follow the instructions on the website to install the development version of Kwant "from source". There are still a couple of features that we want to implement before making the 1.3 release, but I would imagine that this release would be before the end of the month. Joe