! Simple TB sp3 Nearest Neighboard Silicon
! Parameter tooken from Cohen et al.

program TB
implicit none
complex(8) ,DIMENSION(:,:), ALLOCATABLE :: H,S
real(8) kvec(3),step,rwork(30),eigv(8),vl,vu
complex(8) work(30)
integer n,i,lwork,info,l,j,k,il,iu

lwork=30
ALLOCATE(H(8,8),S(8,8))

open(unit=10,file="bands.dat",status="unknown")
rewind(10)

write(*,*) ' Number of points '
read(5,*) n

step = 0.5d0/real(n)

! L point (1/2,1/2,1/2)
kvec(1)=0.5d0
kvec(2)=0.5d0
kvec(3)=0.5d0

do i=1,n

 call fillmatrix(H,kvec)
 call zheev ('N','U',8,H,8,eigv,work,lwork,rwork, info )
 write(10,'(100e15.7)') (eigv(l),l=1,8)
 kvec(1)=kvec(1)-step
 kvec(2)=kvec(2)-step
 kvec(3)=kvec(3)-step
enddo 

! from G point
kvec=0.d0
do i=1,n ! move to X point
 call fillmatrix(H,kvec)
 call zheev ('N','U',8,H,8,eigv,work,lwork,rwork, info )
 write(10,'(100e15.7)') (eigv(l),l=1,8)
 kvec(1)=kvec(1)+2.d0*step
enddo

! Jump to U
kvec(1)=1.d0
kvec(2)=1.d0
kvec(3)=0.d0
!
do i=1,n
 call fillmatrix(H,kvec)
 call zheev ('N','U',8,H,8,eigv,work,lwork,rwork, info )
 write(10,'(100e15.7)') (eigv(l),l=1,8)
 kvec(1) = kvec(1) - 2.d0*step
 kvec(2) = kvec(2) - 2.d0*step
enddo

end

subroutine fillmatrix(H,kvec)   
implicit none
real(8) ess11,epp11,ess12,esp12,esp21,v_pxpx,v_pxpy,epp22,ess22
real(8) sinkx,sinky,sinkz,coskx,cosky,coskz,kx,ky,kz,M_PI
REAL(8),INTENT(IN)    :: kvec(3)
complex(8),INTENT(OUT) :: H(8,8)
complex(8) g0,g1,g2,g3

parameter(M_PI=3.14159265358979323846)


    kx=kvec(1)*M_PI/2.d0
    ky=kvec(2)*M_PI/2.d0
    kz=kvec(3)*M_PI/2.d0

    sinkx=dsin(kx)
    sinky=dsin(ky)
    sinkz=dsin(kz)
    coskx=dcos(kx)
    cosky=dcos(ky)
    coskz=dcos(kz)


    g0 = CMPLX(coskx*cosky*coskz ,-sinkx*sinky*sinkz)
    g1 = CMPLX(-coskx*sinky*sinkz,sinkx*cosky*coskz)
    g2 = CMPLX(-sinkx*cosky*sinkz,coskx*sinky*coskz)
    g3 = CMPLX(-sinkx*sinky*coskz,coskx*cosky*sinkz)
    
    ess11  = 0.d0   
    ess22  = ess11
    epp11  = 7.20 - ess11
    epp22  = epp11
    ess12  = -8.13
    esp12  =  5.88
    esp21  = esp12
    v_pxpx   = 3.17
    v_pxpy   = 7.51  

H=0.d0
    
H(1,1) = CMPLX(ess11) ! S1xS1 
H(2,2) = CMPLX(ess22) ! S2xS2
H(3,3) = CMPLX(epp11) ! Px1 x Px1
H(4,4) = H(3,3)       
H(5,5) = H(3,3)
H(6,6) = CMPLX(epp22) ! Px2 x Px2
H(7,7) = H(6,6)
H(8,8) = H(6,6)

! the off-diagonal parts
H(1,2) = ess12*g0   !s1*s1
H(2,1) = conjg(H(1,2))
! S2 x P1

H(3,2) = -esp12*g1 
H(2,3) = conjg(H(3,2))
H(4,2) = -esp12*g2
H(2,4) = conjg(H(4,2))
H(5,2) = -esp12*g3
H(2,5) = conjg(H(5,2))

! S1 x P2

H(1,6) = esp21*g1
H(6,1) = conjg(H(1,6)) 
H(1,7) = esp21*g2
H(7,1) = conjg(H(1,7)) 
H(1,8) = esp21*g3
H(8,1) = conjg(H(1,8)) 

! P1xP2
H(3,6) = v_pxpx*g0
H(6,3) = conjg(H(3,6))
H(3,7) = v_pxpy*g3
H(7,3) = conjg(H(3,7))
H(3,8) = v_pxpy*g2
H(8,3) = conjg(H(3,8))

H(4,6) = v_pxpy*g3
H(6,4) = conjg(H(4,6))
H(4,7) = v_pxpx*g0
H(7,4) = conjg(H(4,7))
H(4,8) = v_pxpy*g1
H(8,4) = conjg(H(4,8))

H(5,6) = v_pxpy*g2
H(6,5) = conjg(H(5,6))
H(5,7) = v_pxpy*g1
H(7,5) = conjg(H(5,7))
H(5,8) = v_pxpx*g0
H(8,5) = conjg(H(5,8))
end