!*************************************************************************
Program PW
!*************************************************************************
! PW calculates the band structure e(k) for a crystal using a 
! plane wave basis and a fixed non-self-consistent potential. 
! Several choices of empirical pseudopotentials are provided.
! Other choices can be added by modifiying the function "atomPotential"
!    in "atomPotentialMod.90"
! Options:
! Plotting the bands along specified directions in k-space
! Standard diagonalization or conjugate gradient band-by-band minimization
!
! Written by Todd Beaudet of the University of Illinois, based upon programs
! written by W. Mattson, D. Das, K. Glassford, I. Souza, and R. M. Martin
! Contributions also by T. Burns
!
! Note: atomic units (a.u.) are used throughout the program.
!
!------------------------------------------------------------------------

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

  use TagHandlerMod,   only : OpenInputFile, CloseInputFile, TagValue
  use StructureMod,    only : StructureT, StructureInit, StructureDestroy
  use kPointsMod,      only : kPointsT, KPointsInit, KPointsDestroy
  use graphMod,        only : graphT, PlotBands
  use hamiltonianMod,  only : hamiltonianArrayT, HInit, HDestroy, &
                                HDiagonalize
  use eigenStatesMod,  only : eigenStatesT, EigenStatesInit, &
                                EigenStatesDestroy, PrintEigenvalues

  ! PW modules  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  use conjGradMod,     only : conjGradT, conjGradInit, conjGradDestroy, &
                                CGBandByBandMinimize
  use pwHamMod,        only : hamInfoT, HamInfoInit, HamInfoDestroy, &
                                GetBasisAndKE, HGenerate
  use GridMod,         only : GridT, GridInit, GridDestroy, &
                                PutPotentialOnGrid

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

  implicit none

  type(StructureT),        pointer :: structure
  type(kPointsT),          pointer :: kPoints
  type(graphT),            pointer :: graphinfo
  type(hamInfoT),          pointer :: hamInfo
  type(eigenStatesT),      pointer :: eigenStates
  type(hamiltonianArrayT), pointer :: hamiltonian
  type(conjGradT),         pointer :: conjGrad
  type(GridT),             pointer :: grid
  logical                          :: useCGMethod = .false.     ! default
  integer                          :: kpt


  call OpenInputFile()
  call TagValue('DiagonalizationSwitch', logic = useCGMethod)
  if (.not. useCGMethod) print *, 'Using LAPACK for diagonalization'
  if (useCGMethod)  print *, 'Using CG FFT iterative diagonalization'

  call StructureInit(structure)
  call KPointsInit(kPoints, graphInfo, structure%ndim, structure%lattice)
  call HamInfoInit(hamInfo, structure, kPoints)
  call EigenStatesInit(eigenStates, kPoints%numBands, kPoints%numKPoints, &
                       hamInfo%maxBasisSize, -1.0D0, 1)

  if (useCGMethod) then
    call ConjGradInit(conjGrad, hamInfo%maxBasisSize)
    call GridInit(grid, hamInfo)               ! set up grid and maps
    call PutPotentialOnGrid(grid, hamInfo)     ! get real space V on grid
  else
    call HInit(hamiltonian, hamInfo%maxBasisSize)
  end if

  do kpt = 1, kPoints%numKPoints 
    call GetBasisAndKE(hamInfo, eigenStates, kPoints, structure, kpt)
    if (useCGMethod) then
      call CGBandByBandMinimize(conjGrad, grid, hamInfo, eigenStates, kpt)
    else
      call HGenerate(hamInfo, eigenStates, kpt, &
                     hamiltonian%h, hamiltonian%sqH, hamiltonian%hSize)
      call HDiagonalize(hamiltonian, eigenStates, kpt)
    end if
    print *, kpt, 'completed out of a total of ', kPoints%numKPoints
  end do

  if (associated(graphinfo)) &               ! From command line type:
    call PlotBands(graphinfo, eigenStates)   ! gnuplot gnuplot.dat
  if (.not. associated(graphinfo)) &
    call PrintEigenvalues(eigenStates, kPoints%numKPoints)    
  
  ! Clean up memory
  if (useCGMethod) then
    call GridDestroy(grid)
    call ConjGradDestroy(conjGrad)
  else
    call HDestroy(hamiltonian)
  end if
  call EigenStatesDestroy(eigenStates)
  call HamInfoDestroy(hamInfo)
  call KPointsDestroy(kPoints)
  call StructureDestroy(structure)
  call CloseInputFile()

end program PW