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