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) def lead_shape(site): (x, y) = site.pos return (0 <= y < W) syst = kwant.Builder() syst.fill(template, shape, (0, 0)) lead = kwant.Builder(kwant.TranslationalSymmetry([-a, 0])) lead.fill(template, lead_shape, (0, 0)) syst.attach_lead(lead) syst.attach_lead(lead.reversed()) 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() kwant.plotter.bands(syst.leads[0], params=params, 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()