!*************************************************************************
Program TB
!*************************************************************************
! TB calculates the band structure  e(k) for a crystal using a
! empirical tight binding 2-center Slater-Koster form for the matrix elements.
! Parameters can be specified in various ways:
!   A general form for arbitrary angular momenta
!   Harrison's universal matrix elements.
!   Several simple near-neighbor models in 1,2 and 3 d are included;
!   other models can be added by editing one subroutine
! Options:
! Plotting the bands along specified directions in k-space
!
! Written by Nichols Romero and Richard M. Martin of the Univ. of Illinois,
! based upon programs written by William Mattson and Y. Kwon.
! Contributions also by T. Burns and P. Bokes
!
! Note: atomic units (a.u.) are used throughout the program.
!
!------------------------------------------------------------------------

  ! Common modules  - - - - - - - - - - - - - - - - - - - - - - - - - - -

  use sysParams,           only : double, zero
  use TagHandlerMod,       only : OpenInputFile, CloseInputFile
  use structureMod,        only : structureT, structureInit, structureDestroy
  use kPointsMod,          only : kPointsT, kPointsInit, kPointsDestroy
  use graphMod,            only : graphT, PlotBands
  use hamiltonianMod,      only : hamiltonianArrayT, HInit, HDestroy, &
                                  & HDiagonalize,  HPrint
  use eigenStatesMod,      only : eigenStatesT, eigenStatesInit, &
                                  & eigenStatesDestroy, PrintEigenvalues
  use utilsMod

  ! TB modules  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  use tbHamMod,            only : hamInfoT, hamInfoInit, hamInfoDestroy, &
                                  subHamiltonianArrayT, subHInit, &
                                  subHDestroy, subHGenerate, HGenerate

  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  implicit none

  type( StructureT ),           pointer :: structure1
  type( kPointsT ),             pointer :: kPoints
  type( graphT ),               pointer :: graphinfo
  type( hamInfoT ),             pointer :: hamInfo
  type( eigenStatesT ),         pointer :: eigenStates
  type( hamiltonianArrayT ),    pointer :: hamiltonian
  integer         :: kpt
  character(2)    :: outunit

  !*** Specific to the TB method
  type( subHamiltonianArrayT),  pointer :: subHamiltonian

  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  call OpenInputFile()
  call StructureInit(structure1)  !Calls DimensionsInit, LatticeInit, AtomsInit
  call KPointsInit(kPoints, graphInfo, structure1%ndim, structure1%lattice)

  !*** Specific to the TB method
  ! Reads on-site energies, sets up neighbors, orbit info, etc.
  call HamInfoInit(hamInfo, structure1)

  ! Allocate and generate the sub-Hamiltonian array: subH(orbitJ, orbitI, neigh, atomI)
  ! Array of hamiltonian matrix elements in real space
  call SubHInit( subHamiltonian, maxval(hamInfo%orbits), maxval(hamInfo%numChannels), &
               & maxval(hamInfo%neighbors%neighborList(:)%numNeighbors), &
               & structure1%atomBasis%totalNumAtoms )
  call subHGenerate(hamInfo, structure1, subHamiltonian%subH)

  ! Allocate array for H and work arrays of dimension numTotalOrbits
  Print *, 'Allocate array for H and work arrays - dimension = ', &
  & hamInfo%numTotalOrbits
  call HInit( hamiltonian, hamInfo%numTotalOrbits)

  ! Allocate array for eigenStates array
  call EigenStatesInit( eigenStates, kPoints%numBands, kPoints%numKPoints, &
                      & hamiltonian%hDim, -1.0d0, 1)
  eigenStates%eigenVectors = zero

  ! Calculate the Hamiltonian matrix at each k-pt as sum of
  ! phase factors exp(i k xvec) * subH(orbitJ, orbitI, neigh, atomI)
  ! and diagonialize
  do kpt = 1, kPoints%numKPoints
    call HGenerate( hamInfo, structure1, kpt, kPoints, &
                  & subHamiltonian%subH, hamiltonian%h)
    call HPrint( hamiltonian )
    call HDiagonalize( hamiltonian, eigenStates, kpt )
  end do ! kpt

  Print *, 'Eigenvalues calculated and written to file plot.dat'

! FIXME: NO TAG HANDLER DEFINED HERE!
!  call FindTag( tagHandler,'OutputEnergiesInEV', error)
!  if(error == 0) then
!    print *, 'Use OutputEnergiesInEV'
!    outunit = 'eV'
!  else
!    outunit = 'Ha'
!  end if

  ! Create files plot.dat and gnuplot.dat (Display with: gnuplot gnuplot.dat)
  if (associated(graphinfo )) &
    call PlotBands(graphinfo, eigenStates)
  if(.not. associated(graphinfo)) &
    call PrintEigenvalues(eigenStates, 1)

  !  Clean up memory
  call EigenStatesDestroy( eigenStates )
  call SubHDestroy ( subHamiltonian)
  call HDestroy( hamiltonian )
  call HamInfoDestroy( hamInfo )
  call KPointsDestroy( kPoints )
  call StructureDestroy( structure1 )

  call CloseInputFile()

end Program TB