
! Copyright (C) 2018 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 initulr
use modmain
use modulr
use modomp
implicit none
! local variables
integer ik0,ik,ist,jst
integer iq,ig,i,nthd
real(8) tp(2)
! allocatable arrays
integer, allocatable :: idx(:)
real(8), allocatable :: jlgqr(:,:)
! allocate long-range density and magnetisation arrays
if (allocated(rhormt)) deallocate(rhormt)
allocate(rhormt(npcmtmax,natmtot,nqpt))
if (allocated(rhorir)) deallocate(rhorir)
allocate(rhorir(ngtot,nqpt))
if (allocated(magrmt)) deallocate(magrmt)
if (allocated(magrir)) deallocate(magrir)
if (spinpol) then
  allocate(magrmt(npcmtmax,natmtot,ndmag,nqpt))
  allocate(magrir(ngtot,ndmag,nqpt))
end if
if (allocated(rhoqmt)) deallocate(rhoqmt)
allocate(rhoqmt(npcmtmax,natmtot,nqpt))
if (allocated(rhoqir)) deallocate(rhoqir)
allocate(rhoqir(ngtot,nqpt))
if (allocated(magqmt)) deallocate(magqmt)
if (allocated(magqir)) deallocate(magqir)
if (spinpol) then
  allocate(magqmt(npcmtmax,natmtot,ndmag,nqpt))
  allocate(magqir(ngtot,ndmag,nqpt))
  allocate(mommtu(ndmag,natmtot))
  allocate(mommtru(ndmag,natmtot,nqpt))
end if
! allocate potential and magnetic field arrays
if (allocated(vclru)) deallocate(vclru)
allocate(vclru(nqpt))
if (allocated(vsqmt)) deallocate(vsqmt)
allocate(vsqmt(npcmtmax,natmtot,nqpt))
if (allocated(vsqir)) deallocate(vsqir)
allocate(vsqir(ngtot,nqpt))
if (allocated(bfcru)) deallocate(bfcru)
if (allocated(bsqmt)) deallocate(bsqmt)
if (allocated(bsqir)) deallocate(bsqir)
if (spinpol) then
  allocate(bfcru(nqpt,ndmag))
  allocate(bfcmtru(nqpt,natmtot,ndmag))
  allocate(bsqmt(npcmtmax,natmtot,ndmag,nqpt))
  allocate(bsqir(ngtot,ndmag,nqpt))
end if
! G+Q-vector arrays
if (allocated(vgqc)) deallocate(vgqc)
allocate(vgqc(3,ngvec,nqpt))
if (allocated(gqc)) deallocate(gqc)
allocate(gqc(ngvec,nqpt))
if (allocated(ylmgq)) deallocate(ylmgq)
allocate(ylmgq(lmmaxo,ngvec,nqpt))
if (allocated(sfacgq)) deallocate(sfacgq)
allocate(sfacgq(ngvec,natmtot,nqpt))
if (allocated(expqmt)) deallocate(expqmt)
allocate(expqmt(npcmtmax,natmtot,nqpt))
if (allocated(gclgq)) deallocate(gclgq)
allocate(gclgq(ngvec,nqpt))
if (allocated(jlgqrmt)) deallocate(jlgqrmt)
allocate(jlgqrmt(0:lnpsd,ngvec,nspecies,nqpt))
! find the maximum size of the spherical Bessel function array over all species
call findnjcmax
call omp_hold(nqpt,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(jlgqr,ig,tp) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do iq=1,nqpt
  allocate(jlgqr(njcmax,nspecies))
  do ig=1,ngvec
! determine the G+Q-vectors
    vgqc(:,ig,iq)=vgc(:,ig)+vqc(:,iq)
! G+Q-vector length and (theta, phi) coordinates
    call sphcrd(vgqc(:,ig,iq),gqc(ig,iq),tp)
! spherical harmonics for G+Q-vectors
    call genylm(lmaxo,tp,ylmgq(:,ig,iq))
  end do
! generate the spherical Bessel functions j_l(|G+Q|r)
  call genjlgpr(1,gqc(1,iq),jlgqr)
! structure factors for G+Q-vectors
  call gensfacgp(ngvec,vgqc(:,:,iq),ngvec,sfacgq(:,:,iq))
! generate phase factor functions exp(iQ.r) in each muffin-tin
  call genexpmt(1,jlgqr,ylmgq(:,:,iq),ngvec,sfacgq(:,:,iq),expqmt(:,:,iq))
! generate the Coulomb Green's function in G+Q-space
  call gengclgq(.false.,iq,ngvec,gqc(:,iq),gclgq(:,iq))
! compute the spherical Bessel functions j_l(|G+Q|R_mt)
  call genjlgprmt(lnpsd,ngvec,gqc(:,iq),ngvec,jlgqrmt(:,:,:,iq))
  deallocate(jlgqr)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! number of long-range states
nstulr=nstsv*nkpa
! ultracell effective valence charge (this is not necessarily the same as the
! unit cell valence charge times the number of unit cells in the ultracell)
chgvalu=chgval*dble(nkpa)
! allocate eigenvalue array
if (allocated(evalu)) deallocate(evalu)
allocate(evalu(nstulr,nkpt0))
! allocate and initialise occupation number array
if (allocated(occulr)) deallocate(occulr)
allocate(occulr(nstulr,nkpt0))
allocate(idx(nstulr))
do ik0=1,nkpt0
  ik=(ik0-1)*nkpa+1
  call sortidx(nstulr,occsv(1,ik),idx)
  do ist=1,nstulr
    i=idx(nstulr-ist+1)-1
    ik=(ik0-1)*nkpa+i/nstsv+1
    jst=mod(i,nstsv)+1
    occulr(ist,ik0)=occsv(jst,ik)
  end do
end do
deallocate(idx)
return
end subroutine

