
! Copyright (C) 2016 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 genkappa
use modmain
implicit none
! local variables
integer i1,i2,i3,ikpa
integer nkpt0,ik,jk
real(8) v(3)
! allocatable arrays
real(8), allocatable :: vkl0(:,:),vkc0(:,:),wkpt0(:)
if (.not.ultracell) return
! number of kappa-points
nkpa=ngridkpa(1)*ngridkpa(2)*ngridkpa(3)
! allocate the global kappa-point arrays
if (allocated(vkpalu)) deallocate(vkpalu)
allocate(vkpalu(3,nkpa))
if (allocated(vkpal)) deallocate(vkpal)
allocate(vkpal(3,nkpa))
if (allocated(vkpac)) deallocate(vkpac)
allocate(vkpac(3,nkpa))
! generate the kappa-point grids
ikpa=0
do i3=0,ngridkpa(3)-1
  v(3)=dble(i3)/dble(ngridkpa(3))
  do i2=0,ngridkpa(2)-1
    v(2)=dble(i2)/dble(ngridkpa(2))
    do i1=0,ngridkpa(1)-1
      v(1)=dble(i1)/dble(ngridkpa(1))
      ikpa=ikpa+1
! store the kappa-point in ultracell lattice coordinates
      vkpalu(:,ikpa)=v(:)
! kappa-point in Cartesian coordinates
      call r3mv(bvecu,v,vkpac(:,ikpa))
! kappa-point in unit cell lattice coordinates
      call r3mv(binv,vkpac(:,ikpa),vkpal(:,ikpa))
    end do
  end do
end do
! store the existing k-point arrays
allocate(vkl0(3,nkpt),vkc0(3,nkpt),wkpt0(nkpt))
vkl0(:,1:nkpt)=vkl(:,1:nkpt)
vkc0(:,1:nkpt)=vkc(:,1:nkpt)
wkpt0(1:nkpt)=wkpt(1:nkpt)
! number of k+kappa-points
nkpt0=nkpt
nkpt=nkpt0*nkpa
! deallocate and reallocate k-point arrays
deallocate(vkl,vkc,wkpt)
allocate(vkl(3,nkpt),vkc(3,nkpt),wkpt(nkpt))
jk=0
do ik=1,nkpt0
  do ikpa=1,nkpa
    jk=jk+1
    vkl(:,jk)=vkl0(:,ik)+vkpal(:,ikpa)
    vkc(:,jk)=vkc0(:,ik)+vkpac(:,ikpa)
    wkpt(jk)=wkpt0(ik)/dble(nkpa)
  end do
end do
deallocate(vkl0,vkc0,wkpt0)
return
end subroutine

