Dear all Kwant workers!

My problem is about the calculation of energy dispersion of MoS2 in Ref PRB 92, 195402 (2015). I want to use the tight binding Hamiltonian of MoS2 (Eq (2) in Reference)

to obtain energy dispersion in Figure 6(a). Because the honeycomb of monolayer MoS2 is distinct with the graphene, onsite energy of Mo and S is also different. In graphene,

the onsite energy of A and B atoms in unit cell is equivalent and thus the hopping energy of nearest and two nearest atom in different direction is same. However, this is not right

in MoS2 (A----Mo, B----S). So i choose the code to describe onsite energy in MoS2 as follow

def onsite(site):

        return E_M if (site.family == M) else E_X

Is this right?

Furthermore, i solve the energy band in y-direction by considering only one lead (abs(y)<Ly, Ly is the width). Is it correct to use the following code to build the system?

sym = kwant.TranslationalSymmetry(lat.vec((-1, 0)))

sys = kwant.Builder(sym)

My Kwant code cannot obtain the result in Figure 6(a).

Please help me if you are free, Thanks!

My Kwant code:

import kwant

from matplotlib import pyplot

import tinyarray

import math

import numpy

eV, nm=1,1e-9

lamad_M, lamad_X = 0.075*eV, 0.052*eV

e0, e2, ep, ez = -1.094*eV, -1.512*eV, -3.560*eV, -6.886*eV

V_pdd, V_pdp = 3.689*eV, -1.241*eV

V_ddd, V_ddp, V_ddde = -0.895*eV, 0.252*eV, 0.228*eV

V_ppd, V_ppp = 1.225*eV, -0.467*eV;

a0=0.316*nm;

s=1j;

sz=+1;

sqrt_3=math.sqrt(3)

E_M=tinyarray.array([[e0, 0, 0],[0, e2, -s*lamad_M*sz],[0, s*lamad_M*sz, e2]]);

E_X=tinyarray.array([[ep+V_ppp, -s*lamad_X/2*sz, 0],[s*lamad_X/2*sz, ep+V_ppp, 0],[0, 0, ez-V_ppd]]);

t1_MX=1/7*math.sqrt(2/7)*tinyarray.array([[-9*V_pdp+sqrt_3*V_pdd, 3*sqrt_3*V_pdp-V_pdd, 12*V_pdp+sqrt_3*V_pdd],

    [5*sqrt_3*V_pdp+3*V_pdd, 9*V_pdp-sqrt_3*V_pdd, -2*sqrt_3*V_pdp+3*V_pdd],

    [-V_pdp-3*sqrt_3*V_pdd, 5*sqrt_3*V_pdp+3*V_pdd, 6*V_pdp-3*sqrt_3*V_pdd]])

t2_MX=1/7*math.sqrt(2/7)*tinyarray.array([[0, -6*sqrt_3*V_pdp+2*V_pdd, 12*V_pdp+sqrt_3*V_pdd],

    [0, -6*V_pdp-4*sqrt_3*V_pdd, 4*sqrt_3*V_pdp-6*V_pdd],

    [14*V_pdp, 0, 0]])

t3_MX=1/7*math.sqrt(2/7)*tinyarray.array([[9*V_pdp-sqrt_3*V_pdd, 3*sqrt_3*V_pdp-V_pdd, 12*V_pdp+sqrt_3*V_pdd],

    [-5*sqrt_3*V_pdp-3*V_pdd, 9*V_pdp-sqrt_3*V_pdd, -2*sqrt_3*V_pdp+3*V_pdd],

    [-V_pdp-3*sqrt_3*V_pdd, -5*sqrt_3*V_pdp-3*V_pdd, -6*V_pdp+3*sqrt_3*V_pdd]])

t1_MM=1/4*tinyarray.array([[3*V_ddde+V_ddd, sqrt_3/2*(-V_ddde+V_ddd), -3/2*(V_ddde-V_ddd)],

    [sqrt_3/2*(-V_ddde+V_ddd), 1/4*(V_ddde+12*V_ddp+3*V_ddd), sqrt_3/4*(V_ddde-4*V_ddp+3*V_ddd)],

    [-3/2*(V_ddde-V_ddd), sqrt_3/4*(V_ddde-4*V_ddp+3*V_ddd), 1/4*(3*V_ddde+4*V_ddp+9*V_ddd)]])

t2_MM=1/4*tinyarray.array([[3*V_ddde+V_ddd, sqrt_3*(V_ddde-V_ddd), 0],

    [sqrt_3*(V_ddde-V_ddd), V_ddde+3*V_ddd, 0],

    [0                     ,              0,4*V_ddp]])

t3_MM=1/4*tinyarray.array([[3*V_ddde+V_ddd, sqrt_3/2*(-V_ddde+V_ddd), 3/2*(V_ddde-V_ddd)],

    [sqrt_3/2*(-V_ddde+V_ddd), 1/4*(V_ddde+12*V_ddp+3*V_ddd), -sqrt_3/4*(V_ddde-4*V_ddp+3*V_ddd)],

    [3/2*(V_ddde-V_ddd), -sqrt_3/4*(V_ddde-4*V_ddp+3*V_ddd), 1/4*(3*V_ddde+4*V_ddp+9*V_ddd)]])

t1_XX=1/4*tinyarray.array([[3*V_ppp+V_ppd, sqrt_3*(V_ppp-V_ppd), 0],

    [sqrt_3*(V_ppp-V_ppd), V_ppp+3*V_ppd,0],

    [0,                     0,  4*V_ppp]])

t2_XX=tinyarray.array([[V_ppd,  0,  0],

    [0    ,V_ppp, 0],

    [0,  0,   V_ppp]])

t3_XX=1/4*tinyarray.array([[3*V_ppp+V_ppd, -sqrt_3*(V_ppp-V_ppd), 0],

    [-sqrt_3*(V_ppp-V_ppd), V_ppp+3*V_ppd, 0],

    [0,                     0        ,4*V_ppp]])

lat = kwant.lattice.honeycomb()

M, X = lat.sublattices

def make_lead(Lx=1, Ly=48):

    def Rectangle(pos):

        x, y = pos

        return abs(y)<Ly

    def onsite(site):

        return E_M if (site.family == M) else E_X

    #sym = kwant.TranslationalSymmetry((-Lx, 0))

    #sym = kwant.TranslationalSymmetry((-3, 0))

    sym = kwant.TranslationalSymmetry(lat.vec((-1, 0)))

    sys = kwant.Builder(sym)

    #sys = kwant.Builder()

    sys[lat.shape(Rectangle, (0, 0))] = onsite

    #hopping_onsit_M=((0, 0),M,M)

    #hopping_onsit_X=((0, 0),X,X)

    #sys[kwant.builder.HoppingKind(*hopping_onsit_M)] = E_M##########

    #sys[kwant.builder.HoppingKind(*hopping_onsit_X)] = E_X#########

    hopping2_M =((1, 0), M, M)

    hopping1_M =((-1, 1), M, M)

    hopping3_M =((0, -1), M, M)

    sys[kwant.builder.HoppingKind(*hopping1_M)] = t1_MM

    sys[kwant.builder.HoppingKind(*hopping2_M)] = t2_MM

    sys[kwant.builder.HoppingKind(*hopping3_M)] = t3_MM

    hopping2_X =((1, 0), X, X)

    hopping1_X =((-1, 1), X, X)

    hopping3_X =((0, -1), X, X)

    sys[kwant.builder.HoppingKind(*hopping1_X)] = t1_XX

    sys[kwant.builder.HoppingKind(*hopping2_X)] = t2_XX

    sys[kwant.builder.HoppingKind(*hopping3_X)] = t3_XX

    hopping1_MX =((0, 0), M, X)

    hopping2_MX =((0, 1), M, X)

    hopping3_MX =((-1, 1), M, X)

    sys[kwant.builder.HoppingKind(*hopping1_MX)] = t1_MX

    sys[kwant.builder.HoppingKind(*hopping2_MX)] = t2_MX

    sys[kwant.builder.HoppingKind(*hopping3_MX)] = t3_MX

    def family_colors(site):

        return 0 if (site.family == M)  else 1

    def hopping_lw(site1, site2):

        return 0.05 if site1.family == site2.family else 0.1

    def hopping_color(site1,site2):

         return 'g' if site1.family==site2.family else 'b'

    kwant.plot(sys,site_color=family_colors,site_lw=0.05, hop_lw=hopping_lw,hop_color=hopping_color,colorbar=False)

    return sys

def main():

    lead = make_lead().finalized()

    kwant.plotter.bands(lead, show=False)

    pyplot.xlabel("momentum [(lattice constant)^-1]")

    pyplot.ylabel("energy [eV]")

    pyplot.show()

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