!{\src2tex{textfont=tt}}
!!****m* ABINIT/m_paral_pert
!! NAME
!!  m_paral_pert
!!
!! FUNCTION
!!  This module provides routines to manage the parallelisation/distribution
!!  over perturbations
!!
!! COPYRIGHT
!! Copyright (C) 2013-2014 ABINIT group (MT,FJ)
!! 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.
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

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

#include "abi_common.h"

MODULE m_paral_pert

 use defs_basis
 use defs_abitypes
 use m_profiling
 use m_errors

 use m_xmpi,      only : xcomm_rank, xcomm_size, xmpi_allgather
 use m_paw_an,    only : paw_an_type, paw_an_destroy, paw_an_redistribute
 use m_paw_ij,    only : paw_ij_type, paw_ij_destroy, paw_ij_redistribute
 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_destroy, pawfgrtab_redistribute
 use m_pawrhoij,  only : pawrhoij_type, pawrhoij_destroy, pawrhoij_redistribute
 use m_paral_atom,only : get_proc_atmtab

 implicit none

 private

!public procedures.
 public :: set_pert_comm
 public :: unset_pert_comm
 public :: set_pert_paw
 public :: unset_pert_paw

!private procedures.
 private :: get_exchatom_list

!!***

CONTAINS
!!***

!----------------------------------------------------------------------

!!****f* m_paral_pert/set_pert_comm
!! NAME
!! set_pert_comm
!!
!! FUNCTION
!! Set the MPI communicators over the perturbed cell
!!
!! COPYRIGHT
!! Copyright (C) 2013-2014 ABINIT group (FJ,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
!!  nppert=number of processes for the parallelization over perturbations
!!
!! SIDE EFFECTS
!!  mpi_enreg=informations about MPI parallelization
!!
!! PARENTS
!!      loper3
!!
!! CHILDREN
!!      get_proc_atmtab,xmpi_allgather
!!
!! SOURCE

subroutine set_pert_comm(mpi_enreg,nppert)


!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 'set_pert_comm'
!End of the abilint section

 implicit none

!Arguments ---------------------------------------------
!scalars
 integer,intent(in) :: nppert
 type(MPI_type), intent(inout) :: mpi_enreg
!arrays

!Local variables ---------------------------------------
!scalars
!arrays

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

 if (mpi_enreg%paral_pert==0) return

!New communicators for the perturbed cell
 mpi_enreg%comm_cell =mpi_enreg%comm_cell_pert
 mpi_enreg%me_cell   =xcomm_rank(mpi_enreg%comm_cell)
 mpi_enreg%nproc_cell=mpi_enreg%nproc/nppert

!Adjust other communicators
 mpi_enreg%comm_kpt    =mpi_enreg%comm_cell
 mpi_enreg%me_kpt      =mpi_enreg%me_cell
 mpi_enreg%comm_kptband=mpi_enreg%comm_cell
 mpi_enreg%comm_wvl    =mpi_enreg%comm_cell

end  subroutine set_pert_comm
!!***

!----------------------------------------------------------------------

!!****f* m_paral_pert/unset_pert_comm
!! NAME
!! unset_pert_comm
!!
!! FUNCTION
!! Unset the MPI communicators over the perturbed cell; restore the global communicators.
!!
!! COPYRIGHT
!! Copyright (C) 2013-2014 ABINIT group (FJ,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.
!!
!! SIDE EFFECTS
!!  mpi_enreg=informations about MPI parallelization
!!
!! PARENTS
!!      loper3
!!
!! CHILDREN
!!      get_proc_atmtab,xmpi_allgather
!!
!! SOURCE

subroutine unset_pert_comm(mpi_enreg)


!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 'unset_pert_comm'
!End of the abilint section

 implicit none

!Arguments ---------------------------------------------
!scalars
 type(MPI_type), intent(inout) :: mpi_enreg
!arrays

!Local variables ---------------------------------------
!scalars
!arrays

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

 if (mpi_enreg%paral_pert==0) return

!New communicators for the cell
 mpi_enreg%comm_cell =mpi_enreg%comm_world
 mpi_enreg%nproc_cell=mpi_enreg%nproc
 mpi_enreg%me_cell   =mpi_enreg%me

!Adjust other communicators
 mpi_enreg%comm_kpt    =mpi_enreg%comm_cell
 mpi_enreg%me_kpt      =mpi_enreg%me_cell
 mpi_enreg%comm_kptband=mpi_enreg%comm_cell
 mpi_enreg%comm_wvl    =mpi_enreg%comm_cell

end  subroutine unset_pert_comm
!!***

!----------------------------------------------------------------------

!!****f* m_paral_pert/set_pert_paw
!! NAME
!! set_pert_paw
!!
!! FUNCTION
!! Set PAW data for the parallelization over perturbations:
!!  Set MPI communicator over atomic sites
!!  Redistribute PAW on-site data
!!
!! COPYRIGHT
!! Copyright (C) 2013-2014 ABINIT group (MD,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
!!  dtset <type(dataset_type)>=all input variables for this dataset
!!
!! OUTPUT
!!  paw_an_out(my_natom)<type(paw_an_type)>= --optional--
!!              redistributed PAW arrays given on angular mesh
!!              if not present, paw_an is redistributed "in-place"
!!  paw__ij_out(my_natom)<type(paw_ij_type)>= --optional--
!!              redistributed PAW arrays given on (i,j) channels
!!              if not present, paw_ij is redistributed "in-place"
!!  pawfgrtab_out(my_natom)<type(pawfgrtab_type)>= --optional--
!!              redistributed PAW atomic data given on fine grid
!!              if not present, pawfgrtab is redistributed "in-place"
!!  pawrhoij_out(my_natom)<type(pawrhoij_type)>= --optional--
!!              redistributed PAW rhoij occupancies
!!              if not present, pawrhoij is redistributed "in-place"
!!  old_atmtab=save the indexes of the atoms treated by current proc
!!  old_comm_atom=save the identifier of the MPI communicator
!!
!! SIDE EFFECTS
!!  mpi_enreg=informations about MPI parallelization
!!  my_natom=number of atoms treated by current processor
!!  paw_an(my_natom)<type(paw_an_type)>=PAW arrays given on angular mesh
!!  paw_ij(my_natom)<type(paw_ij_type)>=PAW arrays given on (i,j) channels
!!  pawfgrtab(my_natom)<type(pawfgrtab_type)>=PAW atomic data given on fine grid
!!  pawrhoij(my_natom)<type(pawrhoij_type)>=PAW rhoij occupancies
!!
!! PARENTS
!!      loper3
!!
!! CHILDREN
!!      get_proc_atmtab,xmpi_allgather
!!
!! SOURCE

subroutine set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
&                       paw_an,paw_ij,pawfgrtab,pawrhoij,&
&                       paw_an_out,paw_ij_out,pawfgrtab_out,pawrhoij_out)


!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 'set_pert_paw'
 use interfaces_18_timing
 use interfaces_51_manage_mpi
!End of the abilint section

 implicit none

!Arguments ---------------------------------------------
!scalars
 integer,intent(inout) :: my_natom
 integer,intent(inout) :: old_comm_atom
 type(dataset_type), intent(in) :: dtset
 type(MPI_type), intent(inout) :: mpi_enreg
!arrays
 integer,pointer,intent(out) :: old_atmtab(:)
 type(paw_ij_type),allocatable,target,intent(inout) :: paw_ij(:)
 type(paw_ij_type),optional,pointer,intent(out) :: paw_ij_out(:)
 type(paw_an_type),allocatable,target,intent(inout) :: paw_an(:)
 type(paw_an_type),optional,pointer,intent(out) :: paw_an_out(:)
 type(pawfgrtab_type),allocatable,target,intent(inout) :: pawfgrtab(:)
 type(pawfgrtab_type),optional,pointer,intent(out) :: pawfgrtab_out(:)
 type(pawrhoij_type),allocatable,target,intent(inout) :: pawrhoij(:)
 type(pawrhoij_type),optional,pointer,intent(out) :: pawrhoij_out(:)

!Local variables ---------------------------------------
!scalars
!Type of algo used for communications: 1-brute force, 2-asynchronous
 integer,parameter :: algo_option=2
 logical :: paral_atom
!arrays
 integer,allocatable :: SendAtomProc(:), SendAtomList(:)
 integer,allocatable :: RecvAtomProc(:), RecvAtomList(:)
 real(dp) :: tsec(2)

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

 call timab(593,1,tsec)

 paral_atom=(dtset%natom/=my_natom)

!Default value when parallelization over atomic sites is not activated
 if ((mpi_enreg%paral_pert==0).or.(.not.paral_atom).or.(dtset%usepaw==0)) then
   old_comm_atom=mpi_enreg%comm_atom
   old_atmtab=>mpi_enreg%my_atmtab
   if (present(pawrhoij_out))  pawrhoij_out=>pawrhoij
   if (present(paw_ij_out))    paw_ij_out=>paw_ij
   if (present(paw_an_out))    paw_an_out=>paw_an
   if (present(pawfgrtab_out)) pawfgrtab_out=>pawfgrtab
   return
 end if

!Redefine communicator for atoms and store old communicator
 ABI_ALLOCATE(old_atmtab,(my_natom))
 old_comm_atom=mpi_enreg%comm_atom
 old_atmtab(:)=mpi_enreg%my_atmtab(:)
 call initmpi_atom(dtset,mpi_enreg)
 my_natom=mpi_enreg%my_natom

!If "asynchronous algo", determine lists of atoms to be exchanged
!  between me and other processes
 if (algo_option==2) then
   call timab(594,1,tsec)
   call get_exchatom_list(old_comm_atom,mpi_enreg%comm_atom,old_atmtab,&
&       mpi_enreg%my_atmtab,dtset%natom,mpi_enreg%nproc_pert,SendAtomProc,&
&       SendAtomList,RecvAtomProc,RecvAtomList)
   call timab(594,2,tsec)
 end if

!Redistribute PAW on-site data

!pawrhoij datastructure
 call timab(595,1,tsec)
 if (present(pawrhoij_out)) then
   nullify(pawrhoij_out)
   if (algo_option==1) then
     call pawrhoij_redistribute(pawrhoij,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,pawrhoij_out=pawrhoij_out)
   else if (algo_option==2) then
     call pawrhoij_redistribute(pawrhoij,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,pawrhoij_out=pawrhoij_out, &
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 else
   if (algo_option==1) then
     call pawrhoij_redistribute(pawrhoij,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom)
   else if (algo_option==2) then
     call pawrhoij_redistribute(pawrhoij,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 end if
 call timab(595,2,tsec)

!paw_ij datastructure
 call timab(596,1,tsec)
 if (present(paw_ij_out)) then
   nullify(paw_ij_out)
   if (algo_option==1) then
     call paw_ij_redistribute(paw_ij,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,paw_ij_out=paw_ij_out)
   else if (algo_option==2) then
     call paw_ij_redistribute(paw_ij,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,paw_ij_out=paw_ij_out,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 else
   if (algo_option==1) then
     call paw_ij_redistribute(paw_ij,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom)
   else if (algo_option==2) then
     call paw_ij_redistribute(paw_ij,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 end if
 call timab(596,2,tsec)

!paw_an datastructure
 call timab(597,1,tsec)
 if (present(paw_an_out)) then
   nullify(paw_an_out)
   if (algo_option==1) then
     call paw_an_redistribute(paw_an,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,paw_an_out=paw_an_out)
   else if (algo_option==2) then
     call paw_an_redistribute(paw_an,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,paw_an_out=paw_an_out,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 else
   if (algo_option==1) then
     call paw_an_redistribute(paw_an,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom)
   else if (algo_option==2) then
     call paw_an_redistribute(paw_an,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 end if
 call timab(597,2,tsec)

!pawfgrtab datastructure
 call timab(598,1,tsec)
 if (present(pawfgrtab_out)) then
   nullify(pawfgrtab_out)
   if (algo_option==1) then
     call pawfgrtab_redistribute(pawfgrtab,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,pawfgrtab_out=pawfgrtab_out)
   else if (algo_option==2) then
     call pawfgrtab_redistribute(pawfgrtab,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,pawfgrtab_out=pawfgrtab_out,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 else
   if (algo_option==1) then
     call pawfgrtab_redistribute(pawfgrtab,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom)
   else if (algo_option==2) then
     call pawfgrtab_redistribute(pawfgrtab,old_comm_atom,mpi_enreg%comm_atom,&
&         mpi_atmtab_in=old_atmtab,mpi_atmtab_out=mpi_enreg%my_atmtab,&
&         natom=dtset%natom,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 end if
 call timab(598,2,tsec)

 if (algo_option==2) then
   ABI_DEALLOCATE(SendAtomProc)
   ABI_DEALLOCATE(SendAtomList)
   ABI_DEALLOCATE(RecvAtomProc)
   ABI_DEALLOCATE(RecvAtomList)
 end if

 call timab(593,2,tsec)

end  subroutine set_pert_paw
!!***

!----------------------------------------------------------------------

!!****f* m_paral_pert/unset_pert_paw
!! NAME
!! unset_pert_paw
!!
!! FUNCTION
!! Unset PAW data after the parallelization over perturbations:
!!  Restore MPI communicator over atomic sites
!!  Restore PAW on-site data
!!
!! COPYRIGHT
!! Copyright (C) 2013-2014 ABINIT group (MD,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
!!  dtset <type(dataset_type)>=all input variables for this dataset
!!
!! OUTPUT
!!  paw_an_out(my_natom)<type(paw_an_type)>= --optional--
!!              redistributed PAW arrays given on angular mesh
!!              if not present, paw_an is restored "in-place"
!!  paw__ij_out(my_natom)<type(paw_ij_type)>= --optional--
!!              redistributed PAW arrays given on (i,j) channels
!!              if not present, paw_ij is restored "in-place"
!!  pawfgrtab_out(my_natom)<type(pawfgrtab_type)>= --optional--
!!              redistributed PAW atomic data given on fine grid
!!              if not present, pawfgrtab is restored "in-place"
!!  pawrhoij_out(my_natom)<type(pawrhoij_type)>= --optional--
!!              redistributed PAW rhoij occupancies
!!              if not present, pawrhoij is restored "in-place"
!!  old_atmtab=save the indexes of the atoms treated by current proc
!!  old_comm_atom=save the identifier of the MPI communicator
!!
!! SIDE EFFECTS
!!  mpi_enreg=informations about MPI parallelization
!!  my_natom=number of atoms treated by current processor
!!  paw_an(my_natom)<type(paw_an_type)>=PAW arrays given on angular mesh
!!  paw_ij(my_natom)<type(paw_ij_type)>=PAW arrays given on (i,j) channels
!!  pawfgrtab(my_natom)<type(pawfgrtab_type)>=PAW atomic data given on fine grid
!!  pawrhoij(my_natom)<type(pawrhoij_type)>=PAW rhoij occupancies
!!
!! PARENTS
!!      loper3
!!
!! CHILDREN
!!      get_proc_atmtab,xmpi_allgather
!!
!! SOURCE

subroutine unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
&                       paw_an,paw_ij,pawfgrtab,pawrhoij,&
&                       paw_an_out,paw_ij_out,pawfgrtab_out,pawrhoij_out)


!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 'unset_pert_paw'
 use interfaces_18_timing
 use interfaces_51_manage_mpi
!End of the abilint section

 implicit none

!Arguments ---------------------------------------------
!scalars
 integer,intent(inout) :: my_natom
 integer,intent(inout) :: old_comm_atom
 type(dataset_type), intent(in) :: dtset
 type(MPI_type), intent(inout) :: mpi_enreg
!arrays
 integer,pointer,intent(out) :: old_atmtab(:)
 type(paw_ij_type),allocatable,target,intent(inout) :: paw_ij(:)
 type(paw_ij_type),optional,pointer,intent(out) :: paw_ij_out(:)
 type(paw_an_type),allocatable,target,intent(inout) :: paw_an(:)
 type(paw_an_type),optional,pointer,intent(out) :: paw_an_out(:)
 type(pawfgrtab_type),allocatable,target,intent(inout) :: pawfgrtab(:)
 type(pawfgrtab_type),optional,pointer,intent(out) :: pawfgrtab_out(:)
 type(pawrhoij_type),allocatable,target,intent(inout) :: pawrhoij(:)
 type(pawrhoij_type),optional,pointer,intent(out) :: pawrhoij_out(:)

!Local variables ---------------------------------------
!scalars
!Type of algo used for communications: 1-brute force, 2-asynchronous
 integer,parameter :: algo_option=2
 logical :: exchange,paral_atom
!arrays
 integer,allocatable :: SendAtomProc(:), SendAtomList(:)
 integer,allocatable :: RecvAtomProc(:), RecvAtomList(:)
 real(dp) :: tsec(2)

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

!Nothing to do when parallelization over atomic sites is not activated
 if ((mpi_enreg%paral_pert==0).or.(dtset%usepaw==0)) return
 if (.not.associated(old_atmtab)) return
 paral_atom=(dtset%natom/=size(old_atmtab))
 if (.not.paral_atom) return

!If "asynchronous algo", determine lists of atoms to be exchanged
!  between me and other processes
 exchange=.true.
 if (present(pawrhoij_out).and.present(paw_ij_out).and. &
&    present(paw_an_out)  .and.present(pawfgrtab_out)) exchange=.false.
 if (algo_option==2.and.exchange) then
   call timab(594,1,tsec)
   call get_exchatom_list(mpi_enreg%comm_atom,old_comm_atom,mpi_enreg%my_atmtab,&
&       old_atmtab,dtset%natom,mpi_enreg%nproc_pert,SendAtomProc,&
&       SendAtomList,RecvAtomProc,RecvAtomList)
   call timab(594,2,tsec)
 end if

!Redistribute PAW on-site data (if in-place storage)
!or destroy them (out-of-place storage)

!pawrhoij datastructure
 call timab(595,1,tsec)
 if (present(pawrhoij_out)) then
   call pawrhoij_destroy(pawrhoij_out)
   ABI_DATATYPE_DEALLOCATE(pawrhoij_out)
 else
   if (algo_option==1) then
     call pawrhoij_redistribute(pawrhoij,mpi_enreg%comm_atom,old_comm_atom,&
&         mpi_atmtab_in=mpi_enreg%my_atmtab,mpi_atmtab_out=old_atmtab,&
&         natom=dtset%natom)
   else if (algo_option==2) then
     call pawrhoij_redistribute(pawrhoij,mpi_enreg%comm_atom,old_comm_atom,&
&         mpi_atmtab_in=mpi_enreg%my_atmtab,mpi_atmtab_out=old_atmtab,&
&         natom=dtset%natom,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 end if
 call timab(595,2,tsec)

!paw_ij datastructure
 call timab(596,1,tsec)
 if (present(paw_ij_out)) then
   call paw_ij_destroy(paw_ij_out)
   ABI_DATATYPE_DEALLOCATE(paw_ij_out)
 else
   if (algo_option==1) then
     call paw_ij_redistribute(paw_ij,mpi_enreg%comm_atom,old_comm_atom,&
&         mpi_atmtab_in=mpi_enreg%my_atmtab,mpi_atmtab_out=old_atmtab,&
&         natom=dtset%natom)
   else if (algo_option==2) then
     call paw_ij_redistribute(paw_ij,mpi_enreg%comm_atom,old_comm_atom,&
&         mpi_atmtab_in=mpi_enreg%my_atmtab,mpi_atmtab_out=old_atmtab,&
&         natom=dtset%natom,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 end if
 call timab(596,2,tsec)

!paw_an datastructure
 call timab(597,1,tsec)
 if (present(paw_an_out)) then
   call paw_an_destroy(paw_an_out)
   ABI_DATATYPE_DEALLOCATE(paw_an_out)
 else
   if (algo_option==1) then
     call paw_an_redistribute(paw_an,mpi_enreg%comm_atom,old_comm_atom,&
&         mpi_atmtab_in=mpi_enreg%my_atmtab,mpi_atmtab_out=old_atmtab,&
&         natom=dtset%natom)
   else if (algo_option==2) then
     call paw_an_redistribute(paw_an,mpi_enreg%comm_atom,old_comm_atom,&
&         mpi_atmtab_in=mpi_enreg%my_atmtab,mpi_atmtab_out=old_atmtab,&
&         natom=dtset%natom,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 end if
 call timab(597,2,tsec)

!pawfgrtab datastructure
 call timab(598,1,tsec)
 if (present(pawfgrtab_out)) then
   call pawfgrtab_destroy(pawfgrtab_out)
   ABI_DATATYPE_DEALLOCATE(pawfgrtab_out)
 else
   if (algo_option==1) then
     call pawfgrtab_redistribute(pawfgrtab,mpi_enreg%comm_atom,old_comm_atom,&
&         mpi_atmtab_in=mpi_enreg%my_atmtab,mpi_atmtab_out=old_atmtab,&
&         natom=dtset%natom)
   else if (algo_option==2) then
     call pawfgrtab_redistribute(pawfgrtab,mpi_enreg%comm_atom,old_comm_atom,&
&         mpi_atmtab_in=mpi_enreg%my_atmtab,mpi_atmtab_out=old_atmtab,&
&         natom=dtset%natom,&
&         SendAtomproc=SendAtomProc,SendAtomList=SendAtomList,&
&         RecvAtomProc=RecvAtomProc,RecvAtomList=RecvAtomList)
   end if
 end if
 call timab(598,2,tsec)

!Release some memory
 if (algo_option==2.and.exchange) then
   ABI_DEALLOCATE(SendAtomProc)
   ABI_DEALLOCATE(SendAtomList)
   ABI_DEALLOCATE(RecvAtomProc)
   ABI_DEALLOCATE(RecvAtomList)
 end if

!Restore communicator for atoms
 ABI_DEALLOCATE(old_atmtab)
 call initmpi_atom(dtset,mpi_enreg)
 my_natom=mpi_enreg%my_natom

end  subroutine unset_pert_paw
!!***

!----------------------------------------------------------------------

!!****f* m_paral_atom/get_exchatom_list
!! NAME
!! get_exchatom_list
!!
!! FUNCTION
!! This routine determine the list of atoms to be exchanged by current process and other processes.
!! The initial communicator (mpicomm_in) should be a common ancestor of (mpicomm_in and mpicomm_out);
!! It is splitted in "nbgroup" final communicators (mpicomm_out).
!!
!! COPYRIGHT
!! Copyright (C) 2012-2013 ABINIT group (MD)
!! 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
!!  mpicomm_in=input MPI communicator (over atoms) among which atoms are initially distributed
!!  mpicomm_out=MPI communicator (over atoms) among which we want to distribute atoms
!!  my_atmtab_in= Index of atoms initially treated by current proc (communicator mpicomm_in)
!!  my_atmtab_out= Index of atoms finally treated by current proc (communicator mpicomm_out)
!!  natom = total number of atoms
!!  nbgroup = # of mpicomm_out communicators included in mpicomm_in
!!
!! OUTPUT
!! For current proc :
!! RecvAtomProc(:)= rank of processor from which I expect atom (in mpicomm_in)
!! RecvAtomList(:)= indexes of atoms to be received by me
!!   RecvAtomList(irecv) are the atoms I expect from RecvAtomProc(irecv)
!! SendAtomProc(:)= ranks of process destination of atom (in mpicomm_in)
!! SendAtomList(:)= indexes of atoms to be sent by me
!!   SendAtomList(isend) are the atoms sent to SendAtomProc(isend)
!!
!! PARENTS
!!      m_paral_pert
!!
!! CHILDREN
!!      get_proc_atmtab,xmpi_allgather
!!
!! SOURCE


 subroutine get_exchatom_list(mpicomm_in,mpicomm_out,my_atmtab_in,my_atmtab_out,natom,nbgroup, &
&                             SendAtomProc,SendAtomList,RecvAtomProc,RecvAtomList)


!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 'get_exchatom_list'
!End of the abilint section

 implicit none

!Arguments ---------------------------------------------
!scalars
 integer,intent(in) :: mpicomm_in,mpicomm_out,natom,nbgroup
!arrays
 integer,intent(in) :: my_atmtab_in(:),my_atmtab_out(:)
 integer,allocatable,intent(out) :: RecvAtomProc(:),RecvAtomList(:)
 integer,allocatable,intent(out) :: SendAtomList(:),SendAtomProc(:)

!Local variables ---------------------------------------
!scalars
 integer :: iat_in,iat_out,iatom,iatom_tot,ierr,ii,iproc,iproc_in,iproc_out,iproc_outmax,me_in,me_out
 integer :: my_natom_in,my_natom_out,natom_in,natom_out,nbsend,nbrecv,nproc_in,nproc_out
!arrays
 integer :: buf_int(2)
 integer, allocatable :: atmtab_in(:), atmtab_out(:), buf_int_all(:),ranks(:), sizecomm(:)
 integer,allocatable :: RecvAtomProc_(:),SendAtomProc_(:),SendAtomList_(:),RecvAtomList_(:)

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

 nproc_in=xcomm_size(mpicomm_in)
 me_in=xcomm_rank(mpicomm_in)
 my_natom_in=size(my_atmtab_in)
 nproc_out=xcomm_size(mpicomm_out)
 me_out=xcomm_rank(mpicomm_out)
 my_natom_out=size(my_atmtab_out)

 ABI_ALLOCATE(ranks,(1:nproc_in))
 ABI_ALLOCATE(sizecomm,(1:nproc_in))
 buf_int(1)=me_out; buf_int(2)=nproc_out

 ABI_ALLOCATE(buf_int_all,(2*nproc_in))
 call xmpi_allgather(buf_int,2,buf_int_all,mpicomm_in,ierr)
 do ii=1,nproc_in
   ranks(ii)=buf_int_all(2*ii-1)    !me_out
   sizecomm(ii)=buf_int_all(2*ii)   !nproc_out
 enddo
 ABI_DEALLOCATE(buf_int_all)

!Send
 if (my_natom_in > 0) then
   nbsend=0
   ABI_ALLOCATE(SendAtomList_,(my_natom_in*nbgroup)) !at the maximum
   ABI_ALLOCATE(SendAtomProc_,(my_natom_in*nbgroup)) !at the maximum
   iproc_outmax=natom-1 ! maximum procs that can have atoms (using the way atoms are dispatched, see get_my_atmtab)
   do iproc=1,nproc_in !loops on all procs of all comm_outs (those are included in comm_in)
     iproc_out=ranks(iproc)
     if (iproc_out > iproc_outmax) cycle
     nproc_out=sizecomm(iproc)
     iproc_in=iproc - 1 !the numbering of processor begin with 0 not 1
     call get_proc_atmtab(iproc_out,atmtab_out,natom_out,natom,nproc_out) !get atoms tab that iproc_out needs
     do iat_in=1,my_natom_in
       iatom_tot=my_atmtab_in(iat_in)
       do iat_out=1,natom_out
         if (iatom_tot == atmtab_out(iat_out)) then
           !I have iatom_tot, iatom_tot is needed by iproc_out, so send it
           nbsend=nbsend+1
           SendAtomList_(nbsend)=iatom_tot
           SendAtomProc_(nbsend)=iproc_in  ! send to, sorted by growing process
         end if
       enddo
     enddo
     ABI_DEALLOCATE(atmtab_out)
   enddo
   ABI_ALLOCATE(SendAtomList,(nbsend))
   ABI_ALLOCATE(SendAtomProc,(nbsend))
   SendAtomList(1:nbsend)=SendAtomList_(1:nbsend)
   SendAtomProc(1:nbsend)=SendAtomProc_(1:nbsend)
   ABI_DEALLOCATE(SendAtomList_)
   ABI_DEALLOCATE(SendAtomProc_)
 else
   ABI_ALLOCATE(SendAtomList,(0))
   ABI_ALLOCATE(SendAtomProc,(0))
 end if

!recv
 if (my_natom_out > 0) then !no return before because of xmpi_allgather and of the sending operation
   nbrecv=0
   ABI_ALLOCATE(RecvAtomProc_,(my_natom_out))
   ABI_ALLOCATE(RecvAtomList_,(my_natom_out))
   nproc_in=MIN(nproc_in,natom) !limits of procs than can have atoms (use the way the atoms are dispatched)
   do iproc=1,nproc_in ! loops on procs that can have atoms
     iproc_in=iproc -1
     if (iproc_in/=me_in) then
       !atmtab_in(:) list of atoms that iproc_in has (in the  list of total atoms)
       call get_proc_atmtab(iproc_in,atmtab_in,natom_in,natom,nproc_in)
       do iatom=1,my_natom_out
         iatom_tot=my_atmtab_out(iatom) ! I need iatom_tot
         do iat_in =1, natom_in
           if (iatom_tot == atmtab_in(iat_in)) then  ! iproc_in has it, need to receive this atom from iproc_in
             nbrecv=nbrecv+1
             RecvAtomList_(nbrecv)= iatom_tot
             RecvAtomProc_(nbrecv)=iproc_in ! recv from, sorted by growing process
           end if
         end do
       end do
       ABI_DEALLOCATE(atmtab_in)
     end if
   end do
   ABI_ALLOCATE(RecvAtomList,(nbrecv))
   ABI_ALLOCATE(RecvAtomProc,(nbrecv))
   RecvAtomList(1:nbrecv)=RecvAtomList_(1:nbrecv)
   RecvAtomProc(1:nbrecv)=RecvAtomProc_(1:nbrecv)
   ABI_DEALLOCATE(RecvAtomList_)
   ABI_DEALLOCATE(RecvAtomProc_)
 else
   ABI_ALLOCATE(RecvAtomList,(0))
   ABI_ALLOCATE(RecvAtomProc,(0))
 end if

 ABI_DEALLOCATE(ranks)
 ABI_DEALLOCATE(sizecomm)

end subroutine get_exchatom_list

!----------------------------------------------------------------------

END MODULE m_paral_pert
!!***
