!{\src2tex{textfont=tt}}
!!****f* ABINIT/new_integrate_gamma_tr_lova
!!
!! NAME
!! new_integrate_gamma_tr_lova
!!
!! FUNCTION
!! This routine interpolates gkk onto the the coarse k grid, then
!! integrates the TRANSPORT electron phonon coupling matrices
!! over the kpoints on the fermi surface. A dependency on qpoint
!! remains for gamma_qpt_in/out
!! Copied from new_integrate_gamma
!!
!! COPYRIGHT
!! Copyright (C) 2004-2014 ABINIT group (JPC,MJV)
!! 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
!!
!! OUTPUT
!!   elph_tr_ds%gamma_qpt_trout
!!   elph_tr_ds%gamma_qpt_trin
!!
!! NOTES
!!
!! PARENTS
!!      elphon
!!
!! CHILDREN
!!      complete_gamma,destroy_kptrank,ftgam,ftgam_init,mkqptequiv,wrtout
!!      xmpi_sum
!!
!! SOURCE

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

#include "abi_common.h"

subroutine new_integrate_gamma_tr_lova(elph_ds,elph_tr_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_tr_lova'
 use interfaces_14_hidewrite
 use interfaces_77_ddb, except_this_one => new_integrate_gamma_tr_lova
!End of the abilint section

 implicit none

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

!Local variables-------------------------------
!scalars
 integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ierr,iqpt,isppol, iqpt_fullbz
 integer :: itensor, icomp, jcomp
 integer :: fib1, fib2
 integer :: irec, ik_this_proc
 real(dp) :: etain, etaout
 character(len=500) :: message
 type(kptrank_type) :: kptrank_t
!arrays
 real(dp) :: elvelock(3), elvelockpq(3)
 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)

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

 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_tr_ds%gamma_qpt_trin,(2,9,elph_ds%nbranch*elph_ds%nbranch,elph_ds%nsppol,elph_ds%k_phon%nkpt))
 ABI_CHECK_ALLOC('trying to allocate array elph_tr_ds%gamma_qpt_trin ')
 elph_tr_ds%gamma_qpt_trin = zero

 ABI_ALLOCATE(elph_tr_ds%gamma_qpt_trout,(2,9,elph_ds%nbranch*elph_ds%nbranch,elph_ds%nsppol,elph_ds%k_phon%nkpt))
 ABI_CHECK_ALLOC('trying to allocate array elph_tr_ds%gamma_qpt_trout')
 elph_tr_ds%gamma_qpt_trout = zero

!information
 if (elph_ds%gkqwrite == 0) then
   write (message,'(a)')' new_integrate_gamma_tr : keeping gamma matrices in memory'
   call wrtout(std_out,message,'COLL')
 else if (elph_ds%gkqwrite == 1) then
   write (message,'(a)')' new_integrate_gamma_tr : reading gamma matrices from disk'
   call wrtout(std_out,message,'COLL')
 else
   write (message,'(3a,i3)')' new_integrate_gamma_tr : BUG-',ch10,&
&   ' Wrong value for gkqwrite = ',elph_ds%gkqwrite
   MSG_ERROR(message)
 end if

!allocate temp variables
 ABI_ALLOCATE(tmp_gkk,(2,elph_ds%ngkkband**2,elph_ds%nbranch**2,elph_ds%nsppol,elph_ds%nqpt_full))
 ABI_CHECK_ALLOC(' trying to allocate array tmp_gkkout ')

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

 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_fine%nkpt,nrpt))
 ABI_ALLOCATE(sinkr_fine, (elph_ds%k_fine%nkpt,nrpt))
 call ftgam_init(gprim, elph_ds%k_fine%nkpt,nrpt, elph_ds%k_fine%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%nqptirred
     iqpt_fullbz = elph_ds%qirredtofull(iqpt)
     if (elph_ds%gkqwrite == 0) then
       tmp_gkk(:,:,:,:,iqpt_fullbz) = 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)
     end if
   end do


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

!      NB: we are abusing the complete_gamma routine, for a fixed k and completing
!      the elements of the full gkk matrices.
       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,elph_ds%nsppol
       do ib1=1,elph_ds%ngkkband
         fib1=ib1+elph_ds%minFSband-1
         elvelock(:)=elph_tr_ds%el_veloc(ikpt_phon,fib1,:,isppol)

         do ib2=1,elph_ds%ngkkband
           ibeff=ib2+(ib1-1)*elph_ds%ngkkband
           fib2=ib2+elph_ds%minFSband-1
           elvelockpq(:)= elph_tr_ds%el_veloc(ikpt_phonq,fib2,:,isppol)

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

!          MJV 31/03/2009: Note that the following is valid for any geometry, not just cubic!
!          see eq 5 and 6 of prb 36 4103 (Al-Lehaibi et al 1987)
!          see also Allen PRB 17 3725
!          generalization to tensorial quantities is simple, by keeping the directional
!          references of velock and velockpq as indices.
           do icomp = 1, 3
             do jcomp = 1, 3
               itensor = (icomp-1)*3+jcomp
!              FIXME: could use symmetry i <-> j

               etain  = elvelock(icomp)*elvelockpq(jcomp)
               etaout = elvelock(icomp)*elvelock(jcomp)

               elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,iqpt) = &
&               elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,iqpt) + &
&               tmp_gkq(:,:) &
&               *etain &
&               *elph_ds%gkk_intweight(ib1,ikpt_phon,isppol)*elph_ds%gkk_intweight(ib2,ikpt_phonq,isppol)
               
               elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,iqpt) = &
&               elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,iqpt) + &
&               tmp_gkq(:,:) &
&               *etaout &
&               *elph_ds%gkk_intweight(ib1,ikpt_phon,isppol)*elph_ds%gkk_intweight(ib2,ikpt_phonq,isppol)

             end do
           end do
         end do
       end do

     end do ! isppol
   end do ! iq
 end do ! ik

 call xmpi_sum (elph_tr_ds%gamma_qpt_trout, xmpi_world, ierr)
 call xmpi_sum (elph_tr_ds%gamma_qpt_trin, xmpi_world, ierr)
 call destroy_kptrank (kptrank_t)

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

 ABI_DEALLOCATE(coskr_full)
 ABI_DEALLOCATE(sinkr_full)
 ABI_DEALLOCATE(coskr_fine)
 ABI_DEALLOCATE(sinkr_fine)
!
!normalize tensor with 1/sqrt(v_x**2 * v_y**2)
!
!move the veloc into mka2f_tr_lova, where T dependence is dealt with
!This will cause some slight difference to the results
 if (.true.) then
   do isppol=1, elph_ds%nsppol
     do icomp = 1, 3
       do jcomp = 1, 3
         itensor = (icomp-1)*3+jcomp
         if(abs(elph_tr_ds%FSelecveloc_sq(icomp,isppol))>tol14**2 .and. abs(elph_tr_ds%FSelecveloc_sq(jcomp,isppol))>tol14**2)then
           elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,:) = elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,:) / &
&           sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol))
           elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,:) = elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,:) / &
&           sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol))
         else
!          XG120528 Fixed problem with zero velocity
           elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,:)=zero
           elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,:)=zero
         end if
       end do 
     end do 
   end do ! isppol
 end if

!need prefactor of 1/nkpt for each integration over 1 kpoint index.
!NOT INCLUDED IN elph_ds%gkk_intweight
 elph_tr_ds%gamma_qpt_trout = elph_tr_ds%gamma_qpt_trout* elph_ds%occ_factor / elph_ds%k_phon%nkpt
 elph_tr_ds%gamma_qpt_trin  = elph_tr_ds%gamma_qpt_trin * elph_ds%occ_factor / elph_ds%k_phon%nkpt

 write (message,'(2a)')' new_integrate_gamma_tr : transport gamma matrices are calculated ',&
& ' in recip space and for irred qpoints'
 call wrtout(std_out,message,'COLL')

!DEBUG
!write(std_out,*)' new_integrate_gamma_tr_lova: end  elph_tr_ds%gamma_qpt_trin(1,9,1,1,1)=',elph_tr_ds%gamma_qpt_trin(1,9,1,1,1)
!ENDDEBUG


end subroutine new_integrate_gamma_tr_lova
!!***
