Hello,

My name is Enrique and first of all, thanks for the developping of such a useful tool.

My question is related to that one of the edge potentials. I am trying to reproduce the Haldane model with disorder in a rectangular scattering region. My goal is to introduce edge disorder as well as compute the local current density.

I have two problems:

a) Regarding the edge disorder, I cannot remove the disorder in the edges with the leads (I saw your recent previous discussion)

b) Is it possible already to compute directly from Kwant the local current density? If yes, how?

Thanks again,

Enrique


from __future__ import division # so that 1/2 == 0.5, and not 0

from scipy.sparse import spdiags
from math import pi, sqrt, tanh, e
import numpy as np
import kwant
import random
import math 

# For computing eigenvalues
import scipy.sparse.linalg as sla

# For plotting
from matplotlib import pyplot

# Define the graphene lattice
graphene = kwant.lattice.honeycomb()
a, b = graphene.sublattices


edge_potential=10

def make_system(L=10, W=10, pot=0.001):

    t=1
    tt= 0.04 * t
    phi=1*pi/2
    edge=1      # 1 if there is edge disorder
    bulk=0      # 1 if there is bulk disorder
    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
    

        # Modify the scattering region
    if edge==1:
    # change the values of the potential at the edge
     for site in sys.expand(graphene.shape(rec, (0, 1))):
        if sys.degree(site)<=9 and abs(site.pos[0])!=L-4:  #abs(site.pos[0])!=19 execludes the interfaces with the leads
            sys[site]=edge_potential

    elif bulk==1:
    # change the values of the potential at the edge
     for site in sys.expand(graphene.shape(rec, (0, 1))):
        if sys.degree(site)==9 and abs(site.pos[0])!=19:  #abs(site.pos[0])!=19 execludes the interfaces with the leads
            sys[site]=edge_potential


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

    return sys

    
# Call the main function if the script gets executed (as opposed to imported).

if __name__ == '__main__':
    main()t.show()       


def main():

    sys = make_system()
    
    # Check that the system looks as intended.
    
    def family_color(site):
    #print(sys[site])
        if sys[site]==edge_potential:
            return 'green'
        else: return 'black'

    def site_size(site):
        if sys[site]==edge_potential:
            return 0.35
        else:return 0.2
    kwant.plot(sys,site_color=family_color,site_size=site_size)
    # Finalize the system.
    sys = sys.finalized()
    

    # Compute number of neighbors
    i=sys.sites.index(graphene(5,6))
    all_the_neighbors=sys.graph.out_neighbors(i)
    

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