#!/usr/bin/env python

#
# This code is a complete example of the Wiki page
# example.
# However, the way the tight-binding parameters
# are setup is thoroghly explained here.

# Code as shown in the wiki
from __future__ import print_function, division

import sisl
import numpy as np

# Lattice constant
alat = 1.42
length = ( 1.5 ** 2 + 3. / 4. ) ** .5 * alat

# Create carbon atom
C = sisl.Atom(6, R=alat+0.1)

# Create super cell (this is ill-skewed)
sc = sisl.SuperCell([length,length,10,90,90,60], nsc=[3,3,1])
# Rotate super cell
sc = sc.rotate(-30, v=[0,0,1])

# Create graphene unit cell
gr = sisl.Geometry([[0,0,0],[alat,0,0]], atom=C, sc=sc)
# Note this graphene geometry has 2 atoms, A and B lattice atoms.

# Create tight-binding object
tb = sisl.Hamiltonian(gr)

# Denote interactions ranges
#     on-site , nearest neighbour
dR = (0.1     , alat + 0.1)

# We have two atoms in the unit cell
# Lets first create the hopping parameters for
# the first atom

# atom A
# On-site (orthogonal)
tb[0,0] = 0.
# hopping to B (in-cell = 0,0)
tb[0,1] = -2.7
# hopping to B (in neighbour cell = -1,0)
# We first need to get the offset of the atomic
# index for the neighbouring cell
off = gr.sc_index([-1,0,0]) * len(gr)
tb[0,off+1] = -2.7
# The last hopping to B (in neighbour cell = 0, -1)
off = gr.sc_index([0,-1,0]) * len(gr)
tb[0,off+1] = -2.7

# We repeat the same sequence for atom B
tb[1,1] = 0.
# hopping to A (in-cell = 0,0)
tb[1,0] = -2.7
# hopping to A (in neighbour cell = 1,0)
off = gr.sc_index([1,0,0]) * len(gr)
tb[1,off+0] = -2.7
# The last hopping to A (in neighbour cell = 0, 1)
off = gr.sc_index([0,1,0]) * len(gr)
tb[1,off+0] = -2.7

# Now we can calculate the band-structure
nk = 200
kpts = np.linspace(0, 1./3, nk)
eigs = np.empty([nk,2],np.float64)
for ik, k in enumerate(kpts):
    eigs[ik,:] = tb.eigh(k=[k,-k,0])

import matplotlib.pyplot as plt
plt.figure(figsize=(3, 3))
plt.locator_params(nbins=3)
plt.plot(kpts,eigs[:,0])
plt.plot(kpts,eigs[:,1])
plt.show()