[Kwant] tight-binding Hamiltonian

[simon . flury]

Dear kwant users
I was trying to implement the following Hamiltonian: H_j,j+1 = H_j+1,j = -t, 
H_j,j = -0.02*t*cos(k_y + 2pi (phi/W)*j)
for a square lattice of size W*L with attached leads in y-direction. 
Unfortunately it doesn't work. Could anyone give me a hint
how to do this the right way?

Thank you for your answers in advance

Best,

Simon

def make_system(a=1, L=60, W=30):
    offdiagmat = np.zeros((W, W))
    for i in range(W):
        for j in range(W):
            if i == j + 1 or j == i + 1:
                offdiagmat[i][j] = 1


    hamiltonian = "-1.0*offdiagmat - identity(4)*2*0.01*cos(k_y + 
(2*pi*phi/30)*x)"

    hamiltonian = kwant.continuum.sympify(hamiltonian, 
locals=dict(offdiagmat=offdiagmat))

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

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

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

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

    lead = kwant.Builder(kwant.TranslationalSymmetry([0, -a]))
    lead.fill(template, lead_shape, (0, 0))

    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())

    syst = syst.finalized()
    return syst

def analyze_system():
    params = dict(phi=0.2)
    syst = make_system()

    kwant.plotter.bands(syst.leads[0], params=params, momenta=np.linspace(-pi, 
pi, 200), show=False)

    plt.grid(True)
    plt.xlim(-pi, pi)
    plt.xlabel('momentum [1/A]')
    plt.ylabel('energy [eV]')
    plt.show()

def main():
    analyze_system()


if __name__ == '__main__':
    main()