Dear Hu, I looked quickly on your program and it seems in overall fine. Thus, I didn't check the details (the program is heavy to follow) neither checked the details of the article you are mentionning. I wanted just to point out the fact that you are using expressions which may lead you in a mistake concerning the input: math.sqrt(2/7) depending on your python version may give you 0 (2/7==0!) you better use 2./7 for example. check this and see if you get what you want. Hope this helps On Thu, Dec 28, 2017 at 1:59 AM, Gongwei_Hu <gongwei_hu@163.com> wrote:

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

-- Abbout Adel