
! Copyright (C) 2017 T. Mueller, J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine genhmlu(ik0,h)
use modmain
use modulr
use modomp
implicit none
! arguments
integer, intent(in) :: ik0
complex(8), intent(out) :: h(nstulr,nstulr)
! local variables
integer ik,ist,jst,ispn,i,j
integer ikpa,jkpa,iq,ifq
integer iv(3),igk,ifg,nthd
! automatic arrays
integer idx(nstsv)
! allocatable arrays
complex(8), allocatable :: apwalm(:,:,:,:),evecfv(:,:),evecsv(:,:)
complex(8), allocatable :: wfmt(:,:,:,:),wfir(:,:,:),wfgk(:,:,:)
complex(8), allocatable :: vmat(:,:)
allocate(evecfv(nmatmax,nstfv),evecsv(nstsv,nstsv))
! central k-point
ik=(ik0-1)*nkpa+1
! get the ground-state eigenvectors from file for central k-point
call getevecfv('.OUT',ik,vkl(:,ik),vgkl(:,:,:,ik),evecfv)
call getevecsv('.OUT',ik,vkl(:,ik),evecsv)
! find the matching coefficients
allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot))
call match(ngk(1,ik),gkc(:,1,ik),tpgkc(:,:,1,ik),sfacgk(:,:,1,ik),apwalm)
! index to all states
do ist=1,nstsv
  idx(ist)=ist
end do
! calculate the wavefunctions for all states of the central k-point
allocate(wfmt(npcmtmax,natmtot,nspinor,nstsv),wfgk(ngkmax,nspinor,nstsv))
call genwfsv(.false.,.true.,nstsv,idx,ngridg,igfft,ngk(:,ik),igkig(:,:,ik), &
 apwalm,evecfv,evecsv,wfmt,ngkmax,wfgk)
deallocate(apwalm,evecfv,evecsv)
! determine the interstitial wavefunctions in real-space
allocate(wfir(ngtot,nspinor,nstsv))
call omp_hold(nstsv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(ispn,igk,ifg) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do ist=1,nstsv
  do ispn=1,nspinor
    wfir(:,ispn,ist)=0.d0
    do igk=1,ngk(1,ik)
      ifg=igfft(igkig(igk,1,ik))
      wfir(ifg,ispn,ist)=wfgk(igk,ispn,ist)
    end do
    call zfftifc(3,ngridg,1,wfir(:,ispn,ist))
  end do
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! generate the matrix elements for all Q-vectors
call omp_hold(nqpt,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(vmat,ifq,i,j) &
!$OMP PRIVATE(ikpa,jkpa,jst,iv) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do iq=1,nqpt
  allocate(vmat(nstsv,nstsv))
  ifq=iqfft(iq)
  if (spinpol) then
    call genzvbmatk(vsqmt(:,:,ifq),vsqir(:,ifq),bsqmt(:,:,:,ifq), &
     bsqir(:,:,ifq),nstsv,ngk(1,ik),igkig(:,1,ik),wfmt,wfir,wfgk,vmat)
  else
    call genzvmatk(vsqmt(:,:,ifq),vsqir(:,ifq),nstsv,ngk(1,ik),igkig(:,1,ik), &
     wfmt,wfir,wfgk,vmat)
  end if
  j=0
  do jkpa=1,nkpa
    do jst=1,nstsv
      j=j+1
      do ikpa=1,jkpa
        iv(:)=ivq(:,ikpa)-ivq(:,jkpa)
        if (ivqiq(iv(1),iv(2),iv(3)).ne.iq) cycle
        i=(ikpa-1)*nstsv+1
        call zcopy(nstsv,vmat(:,jst),1,h(i,j),1)
      end do
    end do
  end do
  deallocate(vmat)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
deallocate(wfmt,wfir,wfgk)
! add the second-variational eigenvalues to the diagonal
i=0
do ikpa=1,nkpa
  ik=(ik0-1)*nkpa+ikpa
  do ist=1,nstsv
    i=i+1
    h(i,i)=h(i,i)+evalsv(ist,ik)
  end do
end do
return
end subroutine

