Dear all,
I want to discrete the continuous Hamiltonians in PRL 111, 156402 (2013) and calculate the  band structure of the ribbon. So I change the parameters of the code in "Tutorial 2.10. Discretizing continuous Hamiltonians". But I can not obtain the same band structure as shown in Fig.3(b) of the PRL paper. Could you please help me to check my attached code? Thanks!
Yours sincerely
Hosein Khani
###########################
import kwant
import kwant.continuum
import scipy.sparse.linalg
import scipy.linalg
import numpy as np

# For plotting
import matplotlib as mpl
from matplotlib import pyplot as plt

def qsh_system(a=20, L=4000, W=5000):
hamiltonian = """
+ C * identity(4) + M * kron(sigma_0, sigma_z)
- B * (k_x**2 + k_y**2) * kron(sigma_0, sigma_z)
- D * (k_x**2 + k_y**2) * kron(sigma_0, sigma_0)
+ A * k_x * kron(sigma_z, sigma_x)
- A * k_y * kron(sigma_0, sigma_y)
"""

template = kwant.continuum.discretize(hamiltonian, grid=a)

def shape(site):
(x, y) = site.pos
return (0 <= y < W and 0 <= x < L)

(x, y) = site.pos
return (0 <= y < W)

syst = kwant.Builder()
syst.fill(template, shape, (0, 0))

syst = syst.finalized()
return syst
def analyze_qsh():
#######the  follwing  parameters are given to make sure the hamiltonian is the same as
####### PRL 111, 156402 (2013)
AP=0.02851
BP=(0.4381-0.2081)/2.
DP=(0.4381+0.2081)/2.
MP=(-0.19808+0.19153)/2.
CP=(-0.19808-0.19153)/2.

params = dict(A=AP, B=BP, D=DP, M=MP, C=CP)
syst = qsh_system()
momenta=np.linspace(-0.3, 0.3, 201), show=False)

plt.grid()
plt.xlim(-.3, 0.3)
plt.ylim(-0.2, -0.19)
plt.xlabel('momentum [1/A]')
plt.ylabel('energy [eV]')
plt.show()
# get scattering wave functions at E=0
wf = kwant.wave_function(syst, energy=0, params=params)

# prepare density operators
sigma_z = np.array([[1, 0], [0, -1]])
prob_density = kwant.operator.Density(syst, np.eye(4))
spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2)))

# calculate expectation values and plot them
wf_sqr = sum(prob_density(psi) for psi in wf(0))  # states from left lead
rho_sz = sum(spin_density(psi) for psi in wf(0))  # states from left lead

kwant.plotter.map(syst, wf_sqr, show=False)
plt.show()
def main():
analyze_qsh()

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