!{\src2tex{textfont=tt}}
!!****f* ABINIT/new_integrate_gamma
!!
!! NAME
!! new_integrate_gamma
!!
!! FUNCTION
!! This routine interpolates gkk onto the the coarse k grid, then
!! integrates the electron phonon coupling matrix
!! over the kpoints on the fermi surface. A dependency on qpoint
!! remains for gamma_qpt
!!
!! COPYRIGHT
!! Copyright (C) 2004-2014 ABINIT group (BXu)
!! This file is distributed under the terms of the
!! GNU General Public Licence, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! INPUTS
!!   elph_ds = elphon datastructure with data and dimensions
!!      elph_ds%qpt_full = qpoint coordinates
!!      elph_ds%nqptirred = number of irred qpoints
!!      elph_ds%qirredtofull = indexing of the GKK qpoints found
!!
!! OUTPUT
!!   elph_ds = modified elph_ds%gamma_qpt and created elph_ds%gamma_rpt
!!
!! PARENTS
!!      elphon
!!
!! CHILDREN
!!      complete_gamma,destroy_kptrank,ftgam,ftgam_init,mkkptrank,mkqptequiv
!!      wrtout,xmpi_sum
!!
!! SOURCE

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

#include "abi_common.h"

subroutine new_integrate_gamma(elph_ds,Cryst,nrpt,wghatm,rpt,gprim,qpttoqpt)

 use defs_basis
 use defs_elphon
 use m_profiling
 use m_kptrank
 use m_errors
 use m_xmpi

 use m_crystal,   only : crystal_t

!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 'new_integrate_gamma'
 use interfaces_14_hidewrite
 use interfaces_77_ddb, except_this_one => new_integrate_gamma
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 type(elph_type),intent(inout) :: elph_ds
 type(crystal_t),intent(in) :: Cryst
 integer,intent(in) :: nrpt
 integer,intent(in) :: qpttoqpt(2,Cryst%nsym,elph_ds%nqpt_full)
!arrays
 real(dp),intent(in) :: rpt(3,nrpt),wghatm(Cryst%natom,Cryst%natom,nrpt)
 real(dp),intent(in) :: gprim(3,3)

!Local variables-------------------------------
!scalars
 integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ierr,iqpt,isppol!, iqpt_fullbz
 integer :: irec, nbranch,nsppol,ngkkband, ik_this_proc
 character(len=500) :: message
 character(len=fnlen) :: fname
 type(kptrank_type) :: kptrank_t
!arrays
 integer,allocatable :: tmp_FSfullpqtofull(:,:)
 integer,allocatable :: tmp_qpttoqpt(:,:,:)
 real(dp),allocatable :: tmp_gkk(:,:,:,:,:)
 real(dp),allocatable :: tmp_gkk_rpt(:,:,:,:,:)
 real(dp),allocatable :: coskr_full(:,:)
 real(dp),allocatable :: sinkr_full(:,:)
 real(dp),allocatable :: coskr_fine(:,:)
 real(dp),allocatable :: sinkr_fine(:,:)
 real(dp) :: tmp_gkq(2,elph_ds%nbranch*elph_ds%nbranch)

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

 write (message,'(3a)')ch10,' entering new_integrate_gamma ',ch10
 call wrtout(std_out,message,'COLL')

 nsppol   = elph_ds%nsppol
 nbranch  = elph_ds%nbranch
 ngkkband = elph_ds%ngkkband

 ABI_ALLOCATE(tmp_FSfullpqtofull,(elph_ds%k_phon%nkpt,elph_ds%k_phon%nkpt))
 ABI_CHECK_ALLOC("allocating tmp_FSfullpqtofull")

!tmp_qpttoqpt(itim,isym,iqpt) = qpoint index which transforms to iqpt under isym and with time reversal itim.
 ABI_ALLOCATE(tmp_qpttoqpt,(2,Cryst%nsym,elph_ds%k_phon%nkpt))
 ABI_CHECK_ALLOC("allocating tmp_qpttoqpt")

 call mkqptequiv (tmp_FSfullpqtofull,Cryst,elph_ds%k_phon%kpt,elph_ds%k_phon%nkpt,&
& elph_ds%k_phon%nkpt,tmp_qpttoqpt,elph_ds%k_phon%kpt)

 ABI_ALLOCATE(elph_ds%gamma_qpt,(2,nbranch**2,nsppol,elph_ds%k_phon%nkpt))
 ABI_CHECK_ALLOC('out of memory in array elph_ds%gamma_qpt')
 elph_ds%gamma_qpt = zero

 ABI_ALLOCATE(tmp_gkk,(2,ngkkband**2,nbranch**2,nsppol,elph_ds%nqpt_full))
 ABI_CHECK_ALLOC("out of memory in array tmp_gkk")

 ABI_ALLOCATE(tmp_gkk_rpt,(2,ngkkband**2,nbranch**2,nsppol,nrpt))
 ABI_CHECK_ALLOC('out of memory in array tmp_gkk')

 if (elph_ds%gkqwrite == 0) then
   call wrtout(std_out,' new_integrate_gamma : keeping gamma matrices in memory','COLL')
 else if (elph_ds%gkqwrite == 1) then
   fname=trim(elph_ds%elph_base_name) // '_GKKQ'
   write (message,'(2a)')' new_integrate_gamma : reading gamma matrices from file ',trim(fname)
   call wrtout(std_out,message,'COLL')
 else
   write (message,'(a,i0)')' Wrong value for gkqwrite = ',elph_ds%gkqwrite
   MSG_BUG(message)
 end if

 call mkkptrank (elph_ds%k_phon%kpt,elph_ds%k_phon%nkpt,kptrank_t)

 ABI_ALLOCATE(coskr_full, (elph_ds%nqpt_full,nrpt))
 ABI_ALLOCATE(sinkr_full, (elph_ds%nqpt_full,nrpt))
 call ftgam_init(gprim, elph_ds%nqpt_full, nrpt, elph_ds%qpt_full, rpt, coskr_full, sinkr_full)

 ABI_ALLOCATE(coskr_fine, (elph_ds%k_phon%nkpt,nrpt))
 ABI_ALLOCATE(sinkr_fine, (elph_ds%k_phon%nkpt,nrpt))
 call ftgam_init(gprim, elph_ds%k_phon%nkpt, nrpt, elph_ds%k_phon%kpt, rpt, coskr_fine, sinkr_fine)

 do ik_this_proc =1,elph_ds%k_phon%my_nkpt
   ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)

   tmp_gkk = zero
   do iqpt=1,elph_ds%nqpt_full
!    iqpt_fullbz = elph_ds%qirredtofull(iqpt)
     if (elph_ds%gkqwrite == 0) then
       tmp_gkk(:,:,:,:,iqpt) = elph_ds%gkk_qpt(:,:,:,ik_this_proc,:,iqpt)
     else if (elph_ds%gkqwrite == 1) then
       irec = (iqpt-1)*elph_ds%k_phon%my_nkpt+ik_this_proc
       if (ikpt_phon == 1) then
         write (std_out,*) ' new_integrate_gamma  read record ', irec
       end if
!      read (elph_ds%unitgkq,REC=irec) tmp_gkk(:,:,:,:,iqpt_fullbz)
       read (elph_ds%unitgkq,REC=irec) tmp_gkk(:,:,:,:,iqpt)
     end if
   end do


   do ib1=1,ngkkband
     do ib2=1,ngkkband
       ibeff = ib2+(ib1-1)*ngkkband

!      NB: we are abusing the complete_gamma routine, for a fixed k and completing
!      the elements of the full gkk matrices.
!      FIXME complete_gamma may not work here! bxu
       call complete_gamma(Cryst, elph_ds%nbranch, elph_ds%nsppol, elph_ds%nqptirred, elph_ds%nqpt_full,&
&       elph_ds%ep_scalprod, elph_ds%qirredtofull, qpttoqpt, tmp_gkk(:,ibeff,:,:,:))
       do isppol=1,elph_ds%nsppol
         call ftgam(wghatm,tmp_gkk(:,ibeff,:,isppol,:),tmp_gkk_rpt(:,ibeff,:,isppol,:),&
&         Cryst%natom,elph_ds%nqpt_full,nrpt,1,coskr_full, sinkr_full)
       end do
     end do
   end do ! isppol

   do iqpt=1,elph_ds%k_phon%nkpt
     ikpt_phonq = tmp_FSfullpqtofull(ikpt_phon,iqpt)

     do isppol=1,nsppol
       do ib1=1,ngkkband
         do ib2=1,ngkkband
           ibeff = ib2+(ib1-1)*ngkkband

           call ftgam(wghatm,tmp_gkq,tmp_gkk_rpt(:,ibeff,:,isppol,:),&
&           Cryst%natom,1,nrpt,0, coskr_fine(iqpt,:), sinkr_fine(iqpt,:))

           elph_ds%gamma_qpt(:,:,isppol,iqpt) = elph_ds%gamma_qpt(:,:,isppol,iqpt) + &
&           tmp_gkq(:,:)&
&           *elph_ds%gkk_intweight(ib1,ikpt_phon,isppol)*elph_ds%gkk_intweight(ib2,ikpt_phonq,isppol)
!          NOTE: if ngkkband==1 we are using trivial weights since average
!          over bands was done in normsq_gkk (nmsq_gam_sumFS or nmsq_pure_gkk)
         end do ! ib2
       end do ! ib1
     end do ! isppol
   end do ! iqpt
 end do ! ikpt
 ABI_DEALLOCATE(coskr_full)
 ABI_DEALLOCATE(sinkr_full)
 ABI_DEALLOCATE(coskr_fine)
 ABI_DEALLOCATE(sinkr_fine)

 call destroy_kptrank (kptrank_t)
 call xmpi_sum (elph_ds%gamma_qpt, xmpi_world, ierr)

 ABI_DEALLOCATE(tmp_gkk)
 ABI_DEALLOCATE(tmp_gkk_rpt)
 ABI_DEALLOCATE(tmp_FSfullpqtofull)
 ABI_DEALLOCATE(tmp_qpttoqpt)

!need prefactor of 1/nkpt for each integration over 1 kpoint index. NOT INCLUDED IN elph_ds%gkk_intweight
 do iqpt=1,elph_ds%k_phon%nkpt
!  iqpt_fullbz = elph_ds%qirredtofull(iqpt)
!  elph_ds%gamma_qpt(:,:,:,iqpt_fullbz) = elph_ds%gamma_qpt(:,:,:,iqpt_fullbz) / elph_ds%k_phon%nkpt / n0(1) / n0(1)
!  elph_ds%gamma_qpt(:,:,:,iqpt_fullbz) = elph_ds%gamma_qpt(:,:,:,iqpt_fullbz) / elph_ds%k_phon%nkpt / elph_ds%k_phon%nkpt
   elph_ds%gamma_qpt(:,:,:,iqpt) = elph_ds%gamma_qpt(:,:,:,iqpt) * elph_ds%occ_factor / elph_ds%k_phon%nkpt
 end do

 call wrtout(std_out,' new_integrate_gamma: gamma matrices have been calculated for recip space and irred qpoints ',"COLL")

end subroutine new_integrate_gamma
!!***
