!
! Copyright (C) 2016 Quantum ESPRESSO Foundation
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! Author: Mitsuaki Kawamura, U. Tokyo
!----------------------------------------------------------------------------
MODULE elph_scdft_mod
  !--------------------------------------------------------------------------
  !
  IMPLICIT NONE
  !
  PRIVATE
  !
  PUBLIC elph_scdft
  !
  CONTAINS
 !
!----------------------------------------------------------------------------
SUBROUTINE elph_scdft()
  !--------------------------------------------------------------------------
  !
  ! This routine computes the electron-phonon matrix
  ! in the irreducible Brillouin zone and
  ! expand that to whole BZ.
  !
  USE kinds, ONLY : dp
  USE mp_pools, ONLY : npool, inter_pool_comm, my_pool_id
  USE mp_images, ONLY : me_image
  USE io_global,   ONLY : stdout
  USE cell_base, ONLY : at
  USE ions_base, ONLY : nat
  USE start_k, ONLY: nk1, nk2, nk3
  USE qpoint, ONLY : xq, nksq
  USE dynmat, ONLY : dyn, w2
  USE el_phon, ONLY : el_ph_mat, elph_nbnd_min, elph_nbnd_max
  USE control_ph, ONLY : current_iq
  USE modes, ONLY : u
  USE noncollin_module, ONLY : nspin_lsda
  USE io_files, ONLY : prefix, tmp_dir
  !
  INTEGER :: ik, ib, jb, ii, elph_unit, nrcv
  INTEGER :: nbnd_fs
  !
  COMPLEX(dp),ALLOCATABLE :: gep_col(:,:,:), gep(:,:,:,:)
  CHARACTER(100) :: elphname
  INTEGER, EXTERNAL :: find_free_unit
  !
  WRITE(stdout,*) "[elph_scdft]  write elph.dat with symmetries (only on Fermi surfaces)"
  !
  nbnd_fs = elph_nbnd_max - elph_nbnd_min + 1
  !
  WRITE(stdout,*) "[elph_scdft]   Lowest band which contains FS : ", elph_nbnd_min
  WRITE(stdout,*) "[elph_scdft]  Highest band which contains FS : ", elph_nbnd_max
  WRITE(stdout,*) "[elph_scdft]    # of bands which contains FS : ", nbnd_fs
  !
  ! Compute g in each pool
  !
  ALLOCATE(gep(3*nat, elph_nbnd_min:elph_nbnd_max, elph_nbnd_min:elph_nbnd_max, nksq))
  !
  gep(1:3 * nat, elph_nbnd_min:elph_nbnd_max, elph_nbnd_min:elph_nbnd_max, 1:nksq) = 0.0_dp
  !
  DO ik = 1, nksq
     !
     DO ib = elph_nbnd_min, elph_nbnd_max
        DO jb = elph_nbnd_min, elph_nbnd_max
           !
           DO ii = 1, 3*nat
              gep(ii,jb,ib,ik) = DOT_PRODUCT(u(ii,1:3*nat), el_ph_mat(jb,ib,ik,1:3*nat))
           END DO
           gep(1:3*nat,jb,ib,ik) = MATMUL(gep(1:3*nat,jb,ib,ik), dyn(1:3 * nat,1:3 * nat))
           !
        END DO ! jb
     END DO ! ib
     !
  END DO ! ik
  !
  DO ii = 1, 3 * nat
     IF(w2(ii) <= 0_dp) THEN
        gep(ii, elph_nbnd_min:elph_nbnd_max, elph_nbnd_min:elph_nbnd_max, 1:nksq) = 0_dp
     ELSE
        gep(    ii, elph_nbnd_min:elph_nbnd_max, elph_nbnd_min:elph_nbnd_max,1:nksq) &
        & = gep(ii, elph_nbnd_min:elph_nbnd_max, elph_nbnd_min:elph_nbnd_max,1:nksq) &
        &                                  / SQRT(SQRT(w2(ii)) * 2.0_dp)
     END IF
  END DO
  !
  ! Gather El-Ph matrix inter pool
  !
  IF(my_pool_id == 0) THEN
     nrcv = 3 * nat * nbnd_fs * nbnd_fs * nk1*nk2*nk3 * nspin_lsda
     ALLOCATE(gep_col(3*nat*nbnd_fs*nbnd_fs, nk1*nk2*nk3, nspin_lsda))
  ELSE
     nrcv = 1
     ALLOCATE(gep_col(1,1,1))
  END IF
  !
  CALL elph_scdft_gather(gep, 3*nat*nbnd_fs*nbnd_fs*nksq, gep_col, nrcv, &
  &                      my_pool_id, npool, inter_pool_comm)
  !
  ! Write el-ph to file elph.dat
  !
  IF(me_image == 0) THEN
     !
     elph_unit = find_free_unit()
     !
     WRITE(elphname,'(3a,i0)') TRIM(tmp_dir), TRIM(prefix), ".elph", current_iq
     !
     OPEN(elph_unit,file = TRIM(elphname))
     !
     !# of Monkhost-Pack grid
     !
     WRITE(elph_unit,*) nk1, nk2, nk3
     ! 
     !# of band
     !
     WRITE(elph_unit,*) elph_nbnd_min, elph_nbnd_max
     !
     ! q-vector(Crystal cordinate)
     !
     WRITE(elph_unit,*) MATMUL(xq(1:3), at(1:3, 1:3))
     !
     ! # of mode
     !
     WRITE(elph_unit,*) 3 * nat
     !
     ! Frequences[Ryd]
     !
     DO ii = 1, 3 * nat
        WRITE(elph_unit,*) SIGN(SQRT(ABS(w2(ii))), w2(ii))
     END DO
     WRITE(elph_unit,*) ""
     !
     ! Electron-Phonon matrix
     !
     WRITE(elph_unit,'(6e25.15)') &
     &  gep_col(1:3*nat*nbnd_fs*nbnd_fs, 1:nk1*nk2*nk3, 1:nspin_lsda)
     !
     CLOSE(elph_unit)
     !
  END IF ! IF(ionode)
  !
  DEALLOCATE(gep, gep_col)
  !
END SUBROUTINE elph_scdft
!
!------------------------------------------------------------------------
SUBROUTINE elph_scdft_gather(snd,nsnd,rcv,nrcv,mype,npe,comm)
  !----------------------------------------------------------------------
  !
  ! This routine gathers a real matrix to PE 0.
  !
  USE kinds, ONLY : dp
  USE mp, ONLY : mp_sum, mp_gather
  !
  INTEGER,INTENT(IN) :: nsnd, nrcv, mype, npe, comm
  COMPLEX(dp),INTENT(IN) :: snd(nsnd)
  COMPLEX(dp),INTENT(OUT) :: rcv(nrcv)
  !
  INTEGER :: cnt(0:npe - 1), dsp(0:npe - 1), ipe
  !
  cnt(0:npe - 1) = 0
  cnt(mype) = nsnd
  !
  CALL mp_sum(cnt, comm)
  !
  dsp(0) = 0
  DO ipe = 1, npe - 1
     dsp(ipe) = dsp(ipe - 1) + cnt(ipe - 1)
  END DO
  !
  CALL mp_gather(snd(1:nsnd), rcv(1:nrcv), cnt, dsp, 0, comm)
  !
END SUBROUTINE elph_scdft_gather
!
END MODULE elph_scdft_mod
