!{\src2tex{textfont=tt}}
!!****f* ABINIT/posdoppler
!! NAME
!! posdoppler
!!
!! FUNCTION
!! Calculate the momentum distribution annihilating electrons-positron (Doppler broadening)
!!
!! COPYRIGHT
!! Copyright (C) 1998-2014 ABINIT group (JW,GJ,MT)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! INPUTS
!!  atindx1(natom)=index table for atoms, inverse of atindx
!!  cg(2,mcg)=planewave coefficients of wavefunctions.
!!  cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk>
!!                     and each |p_lmn> non-local projector
!!  dimcprj(natom)=array of dimensions of array cprj (not ordered)
!!  dtfil <type(datafiles_type)>=variables related to files
!!   | unpaw=unit number for temporary PAW files
!!  dtset <type(dataset_type)>=all input variables for this dataset
!!   | istwfk=input option parameter that describes the storage of wfs
!!   | mband=maximum number of bands
!!   | mgfft=maximum size of 1D FFTs for the "coarse" grid
!!   | mkmem=number of k points treated by this node.
!!   | mpw=maximum dimensioned size of npw
!!   | natom=number of atoms
!!   | nband=number of bands at each k point
!!   | ngfft=contain all needed information about 3D FFT (coarse grid)
!!   | nkpt=number of k points
!!   | nspden=number of spin-density components
!!   | nspinor=number of spinorial components of the wavefunctions
!!   | nsppol=1 for unpolarized, 2 for spin-polarized
!!   | usepaw=flag for PAW
!!   | use_gpu_cuda=flag for Cuda use
!!   | wtk(=weights associated with various k points
!!  kg(3,mpw*mkmem)=reduced planewave coordinates.
!!  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
!!  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
!!  mpi_enreg= informations about MPI parallelization
!!  npwarr(nkpt)=number of planewaves in basis at this k point
!!  occ(dtset%mband)=occupation numbers
!!  ucvol = unit cell volume
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
!!
!! TODO
!!  nspinor=2 forbidden
!!  NCPP : print a warning because core is missing
!!  force pawusecp=1
!!
!! PARENTS
!!      outscfcv
!!
!! CHILDREN
!!      fourdp,fourwf,pawcprj_alloc,pawcprj_destroy,pawcprj_get
!!      pawcprj_mpi_send,sphereboundary,wrtout,xmpi_recv,xmpi_send
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

subroutine posdoppler(atindx1,cg,cprj,dimcprj,dtfil,dtset,electronpositron,&
&                     kg,mcg,mcprj,mpi_enreg,npwarr,occ,ucvol)

 use m_profiling

 use defs_basis
 use defs_datatypes
 use defs_abitypes
 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_send, pawcprj_destroy
 use m_xmpi
 use m_errors
 use m_electronpositron

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'posdoppler'
 use interfaces_14_hidewrite
 use interfaces_32_util
 use interfaces_52_fft_mpi_noabirule
 use interfaces_53_ffts
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: mcg,mcprj
 real(dp),intent(in) :: ucvol
 type(datafiles_type),intent(in) :: dtfil
 type(dataset_type),intent(in) :: dtset
 type(electronpositron_type),pointer :: electronpositron
 type(MPI_type),intent(inout) :: mpi_enreg
!arrays
 integer,intent(in) :: atindx1(dtset%natom),dimcprj(dtset%natom)
 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
 real(dp),intent(in),target :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
 real(dp),intent(inout),target :: cg(2,mcg)
 type(pawcprj_type),target :: cprj(dtset%natom,mcprj)

!Local variables-------------------------------
!scalars
 integer :: bdtot_index,bdtot_index_pos,cplex,i1,i2,i3,ib,ib_pos
 integer :: ibg,ibg_pos,icg,icg_pos,ierr,ii,ikg,ikg_pos
 integer :: ikpt,ikpt_pos,indx,iorder_cprj,isppol,isppol_pos,istwf_k,istwf_k_pos
 integer :: iwavef,iwavef_pos,ixy,iz,mcg_pos,mcprj_pos,my_nspinor
 integer :: n1,n2,n3,n4,n5,n6,nband_eff_pos,nband_k,nband_k_pos
 integer :: nfft,npw_k,npw_k_pos,nxy
 integer :: option,tag,tim_fourdp=0,tim_fourwf=-36
 logical :: mykpt,mykpt_pos
 real(dp) :: gamma,lambda,units,weight,weight_pos
 character(len=500) :: msg
!arrays
 integer,allocatable :: gbound(:,:),gbound_pos(:,:),kg_k(:,:),kg_k_pos(:,:)
 real(dp),allocatable :: cg_k_pos(:,:),cwaveaug(:,:,:,:),cwaveaug_pos(:,:,:,:)
 real(dp),allocatable :: cwaveg(:,:),cwaveg_pos(:,:),cwaver(:,:),cwaver_pos(:,:)
 real(dp),allocatable :: denpot_dum(:,:,:),fofgout_dum(:,:),occ_k(:),occ_k_pos(:)
 real(dp),allocatable :: rho_contrib(:,:),rho_contrib_g(:,:),rho_moment(:),rho_z(:)
 real(dp),pointer :: cg_pos_ptr(:,:),cg_ptr(:,:),occ_ptr(:),occ_pos_ptr(:)
 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_k_pos(:,:)
 type(pawcprj_type),pointer :: cprj_pos_ptr(:,:),cprj_ptr(:,:)

! *************************************************************************

 DBG_ENTER("COLL")

 if (.not.associated(electronpositron)) then
   msg='  electronpositron variable must be associated !'
   MSG_BUG(msg)
 end if

 if (associated(mpi_enreg%proc_distrb)) then
   do isppol=1,dtset%nsppol
     do ikpt=1,dtset%nkpt
       nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
       if (any(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)/=mpi_enreg%proc_distrb(ikpt,1,isppol))) then
         msg='  proc_distrib cannot be distributed over bands !'
         MSG_BUG(msg)
       end if
     end do
   end do
 end if

!Various initializations
 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
 iorder_cprj=0 ; cplex=2 ; nfft=n1*n2*n3/dtset%ngfft(10)
 ABI_ALLOCATE(rho_contrib,(cplex,nfft))
 ABI_ALLOCATE(rho_contrib_g,(cplex,nfft))
 ABI_ALLOCATE(rho_moment,(nfft))
 rho_moment=zero

!Select electronic and positronic states
 if (electronpositron%particle==EP_ELECTRON) then
   cg_ptr => electronpositron%cg_ep
   cprj_ptr => electronpositron%cprj_ep
   occ_ptr => electronpositron%occ_ep
   cg_pos_ptr => cg
   cprj_pos_ptr => cprj
   occ_pos_ptr => occ
 end if
 if (electronpositron%particle==EP_POSITRON) then
   cg_ptr => cg
   cprj_ptr => cprj
   occ_ptr => occ
   cg_pos_ptr => electronpositron%cg_ep
   cprj_pos_ptr => electronpositron%cprj_ep
   occ_pos_ptr => electronpositron%occ_ep
 end if

!===============================================================================
!================ Loop over positronic states ==================================

!LOOP OVER k POINTS
 ibg_pos=0;icg_pos=0;ikg_pos=0;bdtot_index_pos=0;isppol_pos=1
 do ikpt_pos=1,dtset%nkpt

   npw_k_pos=npwarr(ikpt_pos)
   istwf_k_pos=dtset%istwfk(ikpt_pos)
   nband_k_pos=dtset%nband(ikpt_pos+(isppol_pos-1)*dtset%nkpt)
   ABI_ALLOCATE(occ_k_pos,(nband_k_pos))
   occ_k_pos(:)=occ_pos_ptr(1+bdtot_index_pos:nband_k_pos+bdtot_index_pos)
   nband_eff_pos=1
   do ib_pos=1,nband_k_pos
     if (occ_k_pos(ib_pos)>tol8) nband_eff_pos=ib_pos
   end do

   mcg_pos=npw_k_pos*my_nspinor*nband_eff_pos
   ABI_ALLOCATE(cg_k_pos,(2,mcg_pos))
   if (dtset%usepaw==1) then
     mcprj_pos=my_nspinor*nband_eff_pos
     ABI_DATATYPE_ALLOCATE(cprj_k_pos,(dtset%natom,mcprj_pos))
     call pawcprj_alloc(cprj_k_pos,0,dimcprj)
   end if
   ABI_ALLOCATE(gbound_pos,(2*dtset%mgfft+8,2))
   ABI_ALLOCATE(kg_k_pos,(3,npw_k_pos))
   kg_k_pos(:,1:npw_k_pos)=kg(:,1+ikg_pos:npw_k_pos+ikg_pos)
   call sphereboundary(gbound_pos,istwf_k_pos,kg_k_pos,dtset%mgfft,npw_k_pos)

!  Select k-points for current proc
   mykpt_pos=.true.
   mykpt_pos=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_pos,1,nband_k_pos,isppol_pos,mpi_enreg%me_kpt)) !

!  Exchange data (WF components) between procs
   if (mykpt_pos) then
     cg_k_pos(:,1:mcg_pos)=cg_pos_ptr(:,icg_pos+1:icg_pos+mcg_pos)

     if (dtset%usepaw==1) then
       call pawcprj_get(atindx1,cprj_k_pos,cprj_pos_ptr,dtset%natom,1,ibg_pos,ikpt_pos,iorder_cprj,&
&       isppol_pos,dtset%mband,dtset%mkmem,dtset%natom,nband_eff_pos,nband_k_pos,my_nspinor,&
&       dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
     end if
     do ii=0,mpi_enreg%nproc_kpt-1
       if (ii/=mpi_enreg%me_kpt) then
         tag=ikpt_pos+(isppol_pos-1)*dtset%nkpt+2*dtset%nkpt*ii
         call xmpi_send(cg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
         if (dtset%usepaw==1) then
           call pawcprj_mpi_send(dtset%natom,mcprj_pos,dimcprj,1,cprj_k_pos,ii,mpi_enreg%comm_kpt,ierr)
         end if
       end if
     end do
   else
     ii=0;if (associated(mpi_enreg%proc_distrb)) ii=mpi_enreg%proc_distrb(ikpt_pos,1,isppol_pos)
     tag=ikpt_pos+(isppol_pos-1)*dtset%nkpt+2*dtset%nkpt*mpi_enreg%me_kpt
     call xmpi_recv(cg_k_pos,ii,tag,mpi_enreg%comm_kpt,ierr)
     if (dtset%usepaw==1) then
       call pawcprj_mpi_send(dtset%natom,mcprj_pos,dimcprj,1,cprj_k_pos,ii,mpi_enreg%comm_kpt,ierr)
     end if
   end if

   ABI_ALLOCATE(cwaveg_pos,(2,npw_k_pos*my_nspinor))
   ABI_ALLOCATE(cwaveaug_pos,(2,n4,n5,n6))
   ABI_ALLOCATE(cwaver_pos,(cplex,nfft))

!  ===============================================================================
!  ================ Loop over electronic states ==================================

!  LOOP OVER SPINS
   ibg=0;icg=0;ikg=0;bdtot_index=0
   do isppol=1,dtset%nsppol

!    LOOP OVER k POINTS
     ikg=0
     do ikpt=1,dtset%nkpt

       nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)

!      Select k-points for current proc
       mykpt=.true.
       mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me_kpt))
       if (mykpt) then

!        Allocations depending on k-point
         npw_k=npwarr(ikpt)
         istwf_k=dtset%istwfk(ikpt)
         ABI_ALLOCATE(occ_k,(nband_k))
         occ_k(:)=occ_ptr(1+bdtot_index:nband_k+bdtot_index)

         if (dtset%usepaw==1) then
!          Extract cprj for this k-point according to mkmem
           ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natom,my_nspinor*nband_k))
           call pawcprj_alloc(cprj_k,0,dimcprj)
           call pawcprj_get(atindx1,cprj_k,cprj_ptr,dtset%natom,1,ibg,ikpt,iorder_cprj,isppol,&
&           dtset%mband,dtset%mkmem,dtset%natom,nband_k,nband_k,my_nspinor,dtset%nsppol,dtfil%unpaw,&
&           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
         end if
         ABI_ALLOCATE(gbound,(2*dtset%mgfft+8,2))
         ABI_ALLOCATE(kg_k,(3,npw_k))
         kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
         call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k)

         ABI_ALLOCATE(cwaveg,(2,npw_k*my_nspinor))
         ABI_ALLOCATE(cwaveaug,(2,n4,n5,n6))
         ABI_ALLOCATE(cwaver,(cplex,nfft))

!        =========================================================================
!        Loops on positronic bands

         

         do ib_pos=1,nband_eff_pos
           iwavef_pos=(ib_pos-1)*npw_k_pos*my_nspinor
           if (abs(occ_k_pos(ib_pos))>tol8) then

             cwaveg_pos(:,1:npw_k_pos*my_nspinor)= &
&             cg_k_pos(:,iwavef_pos+1:iwavef_pos+npw_k_pos*my_nspinor)


!            Get positronic wave function in real space
             option=0
             weight_pos=occ_k_pos(ib_pos)*dtset%wtk(ikpt_pos)/ucvol

             call fourwf(1,denpot_dum,cwaveg_pos,fofgout_dum,cwaveaug_pos,&
&             gbound_pos,gbound_pos,istwf_k_pos,kg_k_pos,kg_k_pos,&
&             dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw_k_pos,npw_k_pos,&
&             n4,n5,n6,option,mpi_enreg%paral_kgb,tim_fourwf,weight_pos,weight_pos,&
&             use_gpu_cuda=dtset%use_gpu_cuda)

!            Should use fftpac
!            call fftpac(isppol_pos,mpi_enreg,1,cplex*n1,n2,n3,cplex*n4,n5,n6,&
! &                      dtset%ngfft,cwaver_pos,cwaveaug_pos,10)
!            call fftpac(isppol_pos,mpi_enreg,1,cplex*n1,n2,n3,cplex*n4,n5,n6,&
! &                      dtset%ngfft,cwaver_pos,cwaveaug_pos,11)
             indx=1
             do i3=1,n3
               do i2=1,n2
                 do i1=1,n1
                   cwaver_pos(:,indx)=cwaveaug_pos(:,i1,i2,i3)
                   indx=indx+1
                 end do
               end do
             end do

!            =========================================================================
!            Loops on electronic bands

             do ib=1,nband_k
               iwavef=(ib-1)*npw_k*my_nspinor
               if (abs(occ_k(ib))>tol8) then

                 cwaveg(:,1:npw_k*my_nspinor)= &
&                 cg_ptr(:,icg+iwavef+1:icg+iwavef+npw_k*my_nspinor)

!                Get electronic wave function in real space
                 option=0;tim_fourwf=-36
                 weight=occ(ib+bdtot_index)*dtset%wtk(ikpt)/ucvol

                 call fourwf(1,denpot_dum,cwaveg,fofgout_dum,cwaveaug,&
&                 gbound,gbound,istwf_k,kg_k,kg_k,&
&                 dtset%mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,&
&                 n4,n5,n6,option,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,&
&                 use_gpu_cuda=dtset%use_gpu_cuda)

!                Should use fftpac
!                call fftpac(isppol_pos,mpi_enreg,1,cplex*n1,n2,n3,cplex*n4,n5,n6,&
! &                          dtset%ngfft,cwaver,cwaveaug,10)
!                call fftpac(isppol_pos,mpi_enreg,1,cplex*n1,n2,n3,cplex*n4,n5,n6,&
! &                          dtset%ngfft,cwaver,cwaveaug,11)
                 indx=1
                 do i3=1,n3
                   do i2=1,n2
                     do i1=1,n1
                       cwaver(:,indx)=cwaveaug(:,i1,i2,i3)
                       indx=indx+1
                     end do
                   end do
                 end do

                 
!                =========================================================================
!                Compute plane-wave contribution to momentum distribution

!                * Compute Psi^+(r) * Psi^-(r) * gamma(r) in real space
                 gamma=one


                 do ii=1,nfft
                   rho_contrib(1,ii)=gamma*cwaver_pos(1,ii)*cwaver(1,ii)
                 end do
                 gamma=one
                 if (cplex==2) then
                   do ii=1,nfft
                     rho_contrib(1,ii)=rho_contrib(1,ii) &
&                     -gamma*cwaver_pos(2,ii)*cwaver(2,ii)
                     rho_contrib(2,ii)=gamma*(cwaver_pos(1,ii)*cwaver(2,ii) &
&                     +cwaver_pos(2,ii)*cwaver(1,ii))
                   end do
                 end if

!                FFT of (Psi+.Psi-.gamma) to get Intg[(Psi+.Psi-.gamma).exp(-igr)]
                 call fourdp(cplex,rho_contrib_g,rho_contrib,-1,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%paral_kgb,tim_fourdp)

!                Then take the module square and sum into rho and multipy by scaling factor
                 do ii=1,nfft
                   rho_moment(ii)=rho_moment(ii)+ucvol*ucvol*weight*weight_pos*(rho_contrib_g(1,ii)**2+rho_contrib_g(2,ii)**2)

!                End loops over bands
!                    rho_moment(ii)=rho_moment(ii)+weight*weight_pos*(rho_contrib_g(1,ii)**2+rho_contrib_g(2,ii)**2)
!                    rho_moment(ii)=rho_moment(ii)+ucvol*ucvol*weight*weight_pos*(rho_contrib_g(1,ii)**2)
                 end do
               end if
             end do
           end if
         end do

!        End (if mykpt)
         icg = icg + npw_k*my_nspinor*nband_k
         ibg = ibg +       my_nspinor*nband_k
         ikg = ikg + npw_k
         ABI_DEALLOCATE(cwaveaug)
         ABI_DEALLOCATE(cwaveg)
         ABI_DEALLOCATE(cwaver)
         ABI_DEALLOCATE(gbound)
         ABI_DEALLOCATE(kg_k)
         ABI_DEALLOCATE(occ_k)
         if (dtset%usepaw==1) then
           call pawcprj_destroy(cprj_k)
           ABI_DATATYPE_DEALLOCATE(cprj_k)
         end if
       end if
       bdtot_index=bdtot_index+nband_k

!      End loops over k points and spins (electrons)
!    End loops over k points and spins (electrons)
     end do
   end do

!  End loop over k points (positron)
   if (mykpt_pos) then
     icg_pos = icg_pos + npw_k_pos*my_nspinor*nband_k_pos
     ibg_pos = ibg_pos +           my_nspinor*nband_k_pos
     ikg_pos = ikg_pos + npw_k_pos
   end if
   bdtot_index_pos=bdtot_index_pos+nband_k_pos
   ABI_DEALLOCATE(cwaveaug_pos)
   ABI_DEALLOCATE(cwaveg_pos)
   ABI_DEALLOCATE(cwaver_pos)
   ABI_DEALLOCATE(gbound_pos)
   ABI_DEALLOCATE(kg_k_pos)
   ABI_DEALLOCATE(cg_k_pos)
   ABI_DEALLOCATE(occ_k_pos)
   if (dtset%usepaw==1) then
     call pawcprj_destroy(cprj_k_pos)
     ABI_DATATYPE_DEALLOCATE(cprj_k_pos)
   end if
 end do

!Dont forget, in case of parallelism, to sum over the k communicator

!Integrate rho_moment over p
 units=pi*(one/InvFineStruct)**3/Time_Sec/1.e12_dp/electronpositron%posocc
 lambda=sum(rho_moment(1:nfft))*units/ucvol
 write(msg,'(2a,es22.12)') ch10,'###### lambda=',lambda*1000._dp
 call wrtout(std_out,msg,'COLL')

!Compute rho along z axis
 ABI_ALLOCATE(rho_z,(n3))
 rho_z=zero;nxy=n1*n2
 do iz=1,n3
   do ixy=1,nxy
     rho_z(iz)=rho_z(iz)+units*rho_moment((iz-1)*nxy+ixy)
   end do
 end do
 rho_z=rho_z*((two*pi)**two)/(ucvol**(two/three))
 write(msg,'(2a)') ch10,'rho_z:'
 call wrtout(std_out,msg,'COLL')
 do iz=1,n3
   write(msg,'(es22.15)') rho_z(iz)
   call wrtout(std_out,msg,'COLL')
 end do

 ABI_DEALLOCATE(rho_moment)
 ABI_DEALLOCATE(rho_contrib)
 ABI_DEALLOCATE(rho_contrib_g)
 ABI_DEALLOCATE(rho_z)
 DBG_EXIT("COLL")

end subroutine posdoppler
!!***
