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