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>

[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>

[3]
<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>.
if __name__ == '__main__':
   main()

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