!                                                                            
  ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino  
  !                                                                            
  ! 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 .             
  !                                                                            
  !
  !-----------------------------------------------------------------------
  subroutine loadqmesh
  !-----------------------------------------------------------------------
  !
  !  load fine q mesh and distribute among pools
  !
  !-----------------------------------------------------------------------
#include "f_defs.h"
  !
  USE kinds, only : DP
  USE io_global, ONLY : stdout
  use pwcom
  use phcom
  use epwcom, only : filqf
  use el_phon
#ifdef __PARA
  USE io_global, ONLY : ionode_id
  USE mp_global, ONLY : nproc, my_pool_id, nproc_pool,      &
     intra_image_comm, inter_pool_comm, me_pool, root_pool, &
     mpime, intra_pool_comm, kunit, npool  
  USE mp, ONLY : mp_barrier, mp_bcast
#endif
  !
  implicit none
  !
  integer :: ipol, nsize, rest, iq, nxqftot
  real(kind=DP), parameter :: eps = 1.0e-4
  !
  integer, parameter :: &
    iunqf   = 77       ! unit with fine q point mesh (cartesian coord)
  logical :: exst
  !
  ! work variables 
  !
  real(kind=DP), allocatable :: xqf_ (:,:), wqf_ (:)
  !
#ifdef __PARA
  IF (mpime.eq.ionode_id) then
#endif
    !
    ! read fine mesh
    !
    !
    inquire(file = filqf, exist=exst)
    IF (.not. exst) call errore('loadmesh', 'filqf does not exist',1)
    !
    open ( unit = iunqf, file = filqf, status = 'old', form = 'formatted')
    read(iunqf, *) nxqftot
    WRITE( stdout, '(5x,"Size of q point mesh for interpolation: ",i10)' ) nxqftot
    !
#ifdef __PARA
  ENDIF
  !
  CALL mp_bcast (nxqftot, ionode_id, inter_pool_comm)
  CALL mp_bcast (nxqftot, root_pool, intra_pool_comm)
  !
#endif
  !
  !  temporary storage of k points and weights
  !
  allocate( xqf_ (3, nxqftot), wqf_ (nxqftot) )
  !
#ifdef __PARA
  IF (mpime.eq.ionode_id) then
#endif
    !
    ! read q
    !
    DO iq = 1, nxqftot
      !
      read (iunqf, *) (xqf_ (ipol, iq ), ipol=1,3), wqf_ (iq)
      !
      !  bring the q point to crystal coordinates
      CALL cryst_to_cart ( 1, xqf_ (:,iq), at, -1)
      !
    ENDDO
    !
    close ( iunqf )
    !
    !  check the sum of the kpoint weights: 
    !  must be 2 for spin-unpolarized calculations
    !
    IF (abs(sum( wqf_ )-2.d0).gt.eps) call errore &
         ('ephwann','wrong weights in qpoint grid',1)  
    !
    ! phonons carry no spin...
    !
    wqf_ = wqf_ / float(2)
    !
#ifdef __PARA
  ENDIF
  !
  !  scatter the q points of the fine mesh across the pools
  !
  !  this almost certainly doesn't work
  !  should/could probably be rewritten
  !
  CALL errore('loadmesh', '"Parallel qpoints is untested', 1)
  !
  nxqf = nxqftot / npool 
  rest = nxqftot - nxqf * npool 
  IF ( my_pool_id .lt. rest ) nxqf = nxqf + 1
  !
  nsize = 3
  CALL poolscatter (nsize, nxqftot, xqf_ , nxqf, xqf_ )
  nsize = 1
  CALL poolscatter (nsize, nxqftot, wqf_ , nxqf, wqf_ )
  !
  WRITE( stdout, *) ''
  WRITE( stdout, *) 'Parallelization performed over PHONON wavevectors'
  WRITE( stdout, *) ''
  WRITE( stdout, '(5x,"Max number of q points per pool:",7x,i10)' ) nxqf
  !
#else
  nxqf = nxqftot
#endif
  !
  !  now copy the qpoints and weights into smaller arrays
  !  and deallocate the bigger ones
  !
  allocate ( xqf ( 3, nxqf ), wqf( nxqf ) )
  xqf = xqf_ ( :, 1:nxqf )
  wqf = wqf_ ( 1:nxqf )
  deallocate ( xqf_ , wqf_ )
  !
  end subroutine loadqmesh
  !------------------------------------------------------------
  !
  !-----------------------------------------------------------------------
  subroutine loadkmesh_serial
  !-----------------------------------------------------------------------
  !
  !  load fine k mesh (serial case)
  !
  !-----------------------------------------------------------------------
#include "f_defs.h"
  !
  USE kinds, ONLY : DP
  USE io_global, ONLY : stdout
  use pwcom
  use phcom
  use epwcom, ONLY : filkf, filqf
  use el_phon
#ifdef __PARA
  USE io_global, ONLY : ionode_id
  USE mp_global, ONLY : nproc, my_pool_id, nproc_pool,      &
     intra_image_comm, inter_pool_comm, me_pool, root_pool, &
     mpime, intra_pool_comm, kunit
  USE mp, ONLY : mp_barrier, mp_bcast
#endif
  implicit none
  !
  integer :: ik, ikk, ikq, ipol
  real(kind=DP), parameter :: eps = 1.0e-4
  real(kind=DP) :: xxq(3)
  logical :: exst
  !
  integer, parameter :: &
    iunkkf  = 70       ! unit with fine k point mesh (crys coord)
  !
  !
  ! read fine mesh
  !
  !
  inquire(file = filkf, exist=exst)
  IF (.not. exst) call errore('loadmesh', 'filkf does not exist',1)
  !
  open ( unit = iunkkf, file = filkf, status = 'old', form = 'formatted')
  read(iunkkf, *) nksf
  WRITE( stdout, '(5x,"Size of k point mesh for interpolation: ",i10)' ) nksf 
  !
  IF (lgamma) then
    allocate( xkf (3, nksf), wkf (nksf) )
  ELSE
  allocate( xkf (3, 2*nksf), wkf (2*nksf) )
  ENDIF
  !
  IF (lgamma) then
    !
    DO ik = 1, nksf
      read (iunkkf, *) (xkf (ipol, ik), ipol=1,3), wkf (ik)
    ENDDO
    !
    ! bring the points to crystal coordinates
    !
    CALL cryst_to_cart ( nksf, xkf ,at,-1)
    !
  ELSE
    !
    ! bring q to crystal coordinates
    !
    xxq = xq
    CALL cryst_to_cart (1,xxq,at,-1)
    !
    ! read k and generate k+q in crystal coordinates
    !
    DO ik = 1, nksf
      !
      ikk = 2 * ik - 1
      ikq = ikk + 1
      !
      read (iunkkf, *) (xkf (ipol, ikk ), ipol=1,3), wkf (ikk)
      !
      !  bring the k point to crystal coordinates
      CALL cryst_to_cart ( 1, xkf (:,ikk), at, -1)
      !
      xkf (:, ikq) = xkf (:, ikk) + xxq(:)
      wkf ( ikq ) = 0.d0
      !
    ENDDO
    !
    ! redefine nksf to include the k+q points
    !
    nksf = 2 * nksf
    !
  ENDIF
  close ( iunkkf )
  !
  !  check the sum of the kpoint weights: 
  !  must be 2 for spin-unpolarized calculations
  !
  IF (abs(sum( wkf )-2.d0).gt.eps) call errore &
       ('ephwann','wrong weights in kpoint grid',1)  
  !
  nksqf = nksf / 2 
  !
  WRITE( stdout, '(5x,"Max number of k points per pool:",7x,i10)' ) nksf 
  WRITE( stdout, '(5x,"Max number of k point blocks per pool:",7x,i10)' ) nksqf 
  !
  end subroutine loadkmesh_serial
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! The mesh loading needs to be cleaner.  The above code is a nightmare
! Each load should be succinct and have the option of
! 1) read from file
! 2) generate on a grid (either uniform or in the wedge)
! 3) generate a random number of points
!
! The total end result of the mesh loading should be
!  xkf, wkf, xqf, wqf, nxqf, nkstotf

!-----------------------------------------------------------------------
SUBROUTINE loadkmesh_para
!-----------------------------------------------------------------------
!
!  load fine k mesh and distribute among pools
!
!-----------------------------------------------------------------------
#include "f_defs.h"
#ifdef __PARA
  USE io_global, ONLY : ionode_id
  USE mp_global, ONLY : mpime, inter_pool_comm, &
       my_pool_id, npool
  USE mp,        ONLY : mp_bcast, mp_sum
#endif
  USE kinds,     ONLY : DP
  USE io_global, ONLY : stdout
  USE epwcom,    ONLY : filkf, nkf1, nkf2, nkf3, &
       rand_k, rand_nk, mp_mesh
  USE el_phon,   ONLY : nkstotf, nksf, xkf, wkf, &
       nksqf
  USE cell_base,  ONLY : symm_type
  USE symm_base, ONLY : s, sname, t_rev, time_reversal, &
                        cubicsym, hexsym
  USE pwcom,     ONLY : at, bg, ibrav
  implicit none
  !
  logical :: exst
  integer, PARAMETER :: iunkkf  = 70         
  ! unit with fine k point mesh (crys coord)
  real(kind=DP), ALLOCATABLE :: xkf_(:,:), wkf_(:), xkf_tmp(:,:), &
       wkf_tmp(:)
  integer :: ik, ikk, ikq, lower_bnd, upper_bnd, rest, &
       ipol, i, j, k, nrot, tipo
  

#ifdef __PARA
 IF (mpime .eq. ionode_id) THEN
#endif
    INQUIRE(file = filkf, exist=exst)
    IF (exst) THEN ! load from file
       !
       WRITE (stdout, *) '     Using k-mesh file: ', trim(filkf)
       OPEN ( unit = iunkkf, file = filkf, status = 'old', form = 'formatted')
       READ(iunkkf, *) nkstotf 
       !
       ALLOCATE (xkf_ (3, 2*nkstotf), wkf_(2*nkstotf))
       !
       DO ik = 1, nkstotf
          !
          ikk = 2 * ik - 1
          ikq = ikk + 1
          !
          READ (iunkkf, *) xkf_ (:, ikk ), wkf_ (ikk)
          !
          !  bring the k point to crystal coordinates
          CALL cryst_to_cart ( 1, xkf_ (:,ikk), at, -1)
          !
          xkf_ (:, ikq) = xkf_ (:, ikk) 
          wkf_ ( ikq ) = 0.d0
          !
       ENDDO
       !
       ! redefine nkstotf to include the k+q points
       !
       nkstotf = 2 * nkstotf
       !
    ELSEif ( (nkf1.ne.0) .and. (nkf2.ne.0) .and. (nkf3.ne.0) ) THEN ! generate grid
       IF (mp_mesh) THEN
           ! get size of the mp_mesh in the irr wedge 
          WRITE (stdout, '(a,3i4)') '     Using uniform MP k-mesh: ', nkf1, nkf2, nkf3
          IF (ibrav == 4 .or. ibrav == 5) THEN
             CALL hexsym 
             tipo = 2
          ELSEif (ibrav >= 1 .and. ibrav <= 14) THEN
             CALL cubicsym 
             tipo = 1
          ELSEif (ibrav == 0) THEN
             IF (symm_type == 'cubic') THEN
                CALL cubicsym 
                tipo = 1
             ENDIF
             IF (symm_type == 'hexagonal') THEN
                CALL hexsym 
                tipo = 2
             ENDIF
          ELSE
             CALL errore ('loadmesh', 'wrong ibrav', 1)
          ENDIF
          !                                         
          !
          ALLOCATE ( xkf_ (3, 2*nkf1*nkf2*nkf3), wkf_(2*nkf1*nkf2*nkf3) )
          ! the result of this call is just nkstotf
          CALL kpoint_grid ( nrot, time_reversal, s, t_rev, bg, nkf1*nkf2*nkf3, &
               0,0,0, nkf1,nkf2,nkf3, nkstotf, xkf_, wkf_)
          DEALLOCATE (xkf_, wkf_)
          ALLOCATE ( xkf_ (3, 2*nkstotf), wkf_(2*nkstotf)) 
          ALLOCATE (xkf_tmp (3,nkstotf), wkf_tmp(nkstotf))
          CALL kpoint_grid ( nrot, time_reversal, s, t_rev, bg, nkf1*nkf2*nkf3, &
               0,0,0, nkf1,nkf2,nkf3, nkstotf, xkf_tmp, wkf_tmp)
          !  
          ! assign to k and k+q for xkf and wkf 
          ! 
          DO ik = 1, nkstotf
             ikk = 2 * ik - 1
             ikq = ikk + 1
             xkf_(:,ikk) = xkf_tmp(:,ik)
             xkf_(:,ikq) = xkf_tmp(:,ik)
             wkf_(ikk)   = 2.d0 * wkf_tmp(ik)
             wkf_(ikq)   = 0.d0
          ENDDO
          deALLOCATE (xkf_tmp, wkf_tmp)
          !       
          !       
          CALL cryst_to_cart (2*nkstotf, xkf_, at, -1)
          !
          nkstotf = 2 * nkstotf
          !
       ELSE
          !
          WRITE (stdout, '(a,3i4)') '     Using uniform k-mesh: ', nkf1, nkf2, nkf3
          !
          nkstotf = 2 * nkf1 * nkf2 * nkf3
          ALLOCATE ( xkf_ (3, nkstotf), wkf_(nkstotf) )
          wkf_(:) = 0.d0
          DO ik = 1, nkstotf
             wkf_(2*ik-1) = 2.d0/(float(nkstotf/2))
          ENDDO
          DO i = 1, nkf1
             DO j = 1, nkf2
                DO k = 1, nkf3
                   ik = (i-1)*nkf2*nkf3 + (j-1)*nkf3 + k
                   ikk = 2 * ik - 1
                   ikq = ikk + 1
                   xkf_(1, ikk) = float(i-1)/float(nkf1)
                   xkf_(1, ikq) = xkf_(1, ikk) 
                   xkf_(2, ikk) = float(j-1)/float(nkf2)
                   xkf_(2, ikq) = xkf_(2, ikk) 
                   xkf_(3, ikk) = float(k-1)/float(nkf3)
                   xkf_(3, ikq) = xkf_(3, ikk) 
                ENDDO
             ENDDO
          ENDDO
       !
       ENDIF
       !
    ELSEif (rand_k) THEN  ! random points
       ! random grid
       WRITE (stdout, *) '    Using random k-mesh: ', rand_nk
       !
       nkstotf = rand_nk
       ALLOCATE (xkf_ (3, 2*nkstotf), wkf_(2*nkstotf))
       DO ik = 1, nkstotf
          !
          ikk = 2 * ik - 1
          ikq = ikk + 1
          !
          wkf_(ikk) = 2.d0/ float(nkstotf)
          wkf_(ikq) = 0.d0
          CALL random_number(xkf_(:,ikk))
          xkf_(:,ikq) = xkf_(:,ikk)
       ENDDO
       !
       ! redefine nkstotf to include the k+q points
       !
       nkstotf = 2 * nkstotf
       !
    ELSE ! don't know how to get grid
       CALL errore('loadkmesh_para', "Cannot load fine k points", 1)
    ENDIF
#ifdef __PARA
 ENDIF
 CALL mp_bcast (nkstotf, ionode_id, inter_pool_comm)
 !
 !  scatter the k points of the fine mesh across the pools
 !
 nksf = 2 * ( nkstotf / 2 / npool )
 rest = ( nkstotf - nksf * npool ) / 2
 IF (my_pool_id < rest ) THEN
    nksf=nksf+2
    lower_bnd = my_pool_id*nksf + 1
    upper_bnd = lower_bnd + nksf - 1
 ELSE
    lower_bnd = rest*(nksf+2)+(my_pool_id-rest)*nksf + 1
    upper_bnd = lower_bnd + nksf - 1
 ENDIF
 !
 nksqf = nksf / 2 
 IF (.not.ALLOCATEd(xkf_)) ALLOCATE (xkf_(3,nkstotf))
 IF (.not.ALLOCATEd(wkf_)) ALLOCATE (wkf_(  nkstotf))
 CALL mp_bcast(xkf_, ionode_id)
 CALL mp_bcast(wkf_, ionode_id)
 !
#else
 !
 ! In serial the definitions are much easier 
 !
 nksf = nkstotf
 nksqf = nksf / 2 
 lower_bnd = 1
 upper_bnd = nksf
 !
#endif
 !
 !  Assign the weights and vectors to the correct bounds
 !
 ALLOCATE(xkf(3,nksf))
 ALLOCATE(wkf(  nksf))
 xkf(:,:) = xkf_ (:, lower_bnd:upper_bnd)
 wkf(  :) = wkf_ (   lower_bnd:upper_bnd)
 !
 WRITE( stdout, '(5x,"Size of k point mesh for interpolation: ",i10)' ) nkstotf 
 WRITE( stdout, '(5x,"Max number of k points per pool:",7x,i10)' ) nksf 
 !
 IF (ALLOCATED(xkf_)) DEALLOCATE(xkf_)
 IF (ALLOCATED(wkf_)) DEALLOCATE(wkf_)
 !
END SUBROUTINE loadkmesh_para
!-----------------------------------------------------------------------

!SUBROUTINE loadkmesh_serial
!end SUBROUTINE loadkmesh_serial



!-----------------------------------------------------------------------
SUBROUTINE loadqmesh_para
!-----------------------------------------------------------------------
!
!  load fine q mesh and distribute among pools
!
!-----------------------------------------------------------------------
#include "f_defs.h"
#ifdef __PARA
  USE io_global, ONLY : ionode_id
  USE mp_global, ONLY : mpime
#endif
  USE kinds,     ONLY : DP
  USE io_global, ONLY : stdout
  USE epwcom,    ONLY : filqf, nqf1, nqf2, nqf3, rand_q
  USE el_phon,   ONLY : nxqf
  implicit none
  !
  logical :: exst
  integer, PARAMETER :: iunqf  = 77         
  ! unit with fine q point mesh (crys coord)

#ifdef __PARA
 IF (mpime .eq. ionode_id) THEN
#endif

 IF (filqf .ne. '') THEN ! load from file
    inquire(file = filqf, exist=exst)
    IF (.not. exst) CALL errore('loadmesh', 'filqf does not exist',1)
    !
    WRITE (stdout, *) '     Using q-mesh file: ', trim(filqf)
    OPEN ( unit = iunqf, file = filqf, status = 'old', form = 'formatted')
    READ(iunqf, *) nxqf 
    !
 ELSEif ( (nqf1.ne.0) .and. (nqf2.ne.0) .and. (nqf3.ne.0) ) THEN ! generate grid

 ELSEif (rand_q) THEN  ! random points
    
 ELSE ! don't know how to get grid
    CALL errore('loadqmesh_para', "Cannot load fine q points", 1)
 ENDIF

#ifdef __PARA
 ENDIF
#endif
END SUBROUTINE loadqmesh_para
!-----------------------------------------------------------------------



!-----------------------------------------------------------------------
SUBROUTINE loadqmesh_serial
!-----------------------------------------------------------------------
!
!  load fine q mesh on each pool
!
!-----------------------------------------------------------------------
#include "f_defs.h"
  USE kinds, ONLY : DP
  USE io_global, ONLY : stdout
  USE epwcom,    ONLY : filqf, nqf1, nqf2, nqf3, &
       rand_q, rand_nq, mp_mesh
  USE el_phon,   ONLY : xqf, wqf, nksqf, nxqf
  USE cell_base, ONLY : symm_type
  use symm_base, ONLY : s, sname, t_rev, time_reversal, &
       cubicsym, hexsym
  USE pwcom,     ONLY : at, bg, ibrav
  implicit none
  !
  logical :: exst
  integer, PARAMETER :: iunqf  = 77         
  ! unit with fine k point mesh (crys coord)
  integer :: iq , i, j, k, nrot, tipo


  inquire(file = filqf, exist=exst)
  IF (exst) THEN ! load from file
     !
     ! Each pool gets its own copy from the action=read statement
     !
     WRITE (stdout, *) '     Using q-mesh file: ', trim(filqf)
     OPEN ( unit = iunqf, file = filqf, status = 'old', form = 'formatted', action='read')
     READ(iunqf, *) nxqf
     ALLOCATE (xqf(3, nxqf), wqf(nxqf))
     DO iq = 1, nxqf
        READ (iunqf, *) xqf (:, iq), wqf(iq)
     ENDDO
     CLOSE(iunqf)
     !
     ! bring xqf in crystal coordinates
     !
     CALL cryst_to_cart (nxqf, xqf, at, -1)
     !
 ELSEif ( (nqf1.ne.0) .and. (nqf2.ne.0) .and. (nqf3.ne.0) ) THEN ! generate grid
    IF (mp_mesh) THEN
       ! get size of the mp_mesh in the irr wedge 
       WRITE (stdout, '(a,3i4)') '     Using uniform MP q-mesh: ', nqf1, nqf2, nqf3
       IF (ibrav == 4 .or. ibrav == 5) THEN
          CALL hexsym 
          tipo = 2
       ELSEif (ibrav >= 1 .and. ibrav <= 14) THEN
          CALL cubicsym
          tipo = 1
       ELSEif (ibrav == 0) THEN
          IF (symm_type == 'cubic') THEN
             CALL cubicsym
             tipo = 1
          ENDIF
          IF (symm_type == 'hexagonal') THEN
             CALL hexsym
             tipo = 2
          ENDIF
       ELSE
          CALL errore ('loadmesh', 'wrong ibrav', 1)
       ENDIF
       !                                         
       !
       ALLOCATE ( xqf (3, nqf1*nqf2*nqf3), wqf(nqf1*nqf2*nqf3) )
       ! the result of this call is just nkstotf
       CALL kpoint_grid ( nrot, time_reversal, s, t_rev, bg, nqf1*nqf2*nqf3, &
            0,0,0, nqf1,nqf2,nqf3, nxqf, xqf, wqf)
       deALLOCATE ( xqf, wqf) 
       ALLOCATE ( xqf(3, nxqf), wqf(nxqf)) 
       CALL kpoint_grid ( nrot, time_reversal, s, t_rev, bg, nqf1*nqf2*nqf3, &
            0,0,0, nqf1,nqf2,nqf3, nxqf, xqf, wqf)
       wqf = 2 * wqf
       !  
       !       
       CALL cryst_to_cart (nxqf, xqf, at, -1)
       !
    ELSE
       ! currently no offset.  
       ! q's are in crystal coordinates in xqf
       WRITE (stdout, '(a,3i4)') '     Using uniform q-mesh: ', nqf1, nqf2, nqf3
       !
       nxqf = nqf1 * nqf2 * nqf3
       ALLOCATE ( xqf (3, nxqf), wqf(nxqf) )
       wqf(:) = 2.d0/(float(nxqf))
       DO i = 1, nqf1
          DO j = 1, nqf2
             DO k = 1, nqf3
                iq = (i-1)*nqf2*nqf3 + (j-1)*nqf3 + k
                xqf(1, iq) = float(i-1)/float(nqf1)
                xqf(2, iq) = float(j-1)/float(nqf2)
                xqf(3, iq) = float(k-1)/float(nqf3)
             ENDDO
          ENDDO
       ENDDO
       !
    ENDIF
 ELSEif (rand_q) THEN  ! random points
    ! random grid
    WRITE (stdout, *) '    Using random q-mesh: ', rand_nq
    !
    nxqf = rand_nq
    ALLOCATE (xqf(3, nxqf), wqf(nxqf))
    DO iq = 1, nxqf
       !
       wqf(iq) = 2.d0/ float(nxqf)
       CALL random_number(xqf(:,iq))
    ENDDO
    !
 ELSE ! don't know how to get grid
    CALL errore('loadqmesh_serial', "Cannot load fine q points", 1)
 ENDIF
 !
 IF (abs(sum (wqf) - 2.d0) .gt. 1.d-4 ) &
      CALL errore('loadqmesh_serial',"q-point wieghts do not add up to 2", 1)
 !
 WRITE( stdout, '(5x,"Size of q point mesh for interpolation: ",i10)' ) nxqf
 !
end SUBROUTINE loadqmesh_serial
!-----------------------------------------------------------------------