!{\src2tex{textfont=tt}}
!!****m* ABINIT/m_melemts
!! NAME
!!  m_melemts
!!
!! FUNCTION
!!  This module defines an object used as database to store matrix 
!!  elements of several potentials and operators between two Bloch states.
!!  These values are used in the GW part of abinit to evaluate QP energies
!!  using the perturbative approach.
!!
!! COPYRIGHT
!! Copyright (C) 2008-2014 ABINIT group (MG)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!!
!! NOTES
!!  * This module is supposed to be used only in the GW part to facilitate
!!    further developments. The object might change in the future. Thus
!!    contact Matteo Giantomassi if you wish to use this piece of code
!!    for your developments. In particular we might decide to switch 
!!    to ragged arrays Mels(nkibz,nsppol*nspinor**2)%data
!!
!!  * Routines tagged with "@type_name" are tightly connected to the definition of the data type. 
!!    Tightly connected means that the proper functioning of the implementation relies on the 
!!    assumption that the tagged procedure is consistent with the type declaration.
!!    Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure 
!!    that all the tightly connected routines are changed accordingly to accommodate the modification of the data type. 
!!    Typical examples of tightly connected routines are creation, destruction or reset methods.
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

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

#include "abi_common.h"

MODULE m_melemts

 use defs_basis
 use m_errors
 use m_xmpi
 use m_profiling

 use m_fstrings,       only : tolower
 use m_numeric_tools,  only : print_arr

 implicit none
 
 private
 
 public :: reset_mflags
 public :: copy_mflags
 public :: destroy_melements
 public :: init_melements
 public :: herm_melements
 public :: xsum_melements
 public :: print_melements
 public :: zero_melements
!!***

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

!!****t* m_melemts/melements_flags_type
!! NAME
!!
!! FUNCTION
!!  Container for the flags defining the status of the corresponding 
!!  pointer defined in the type melements_type. Possible values are:
!!    * 0 if the correspondending array is not allocated.
!!    * 1 if allocated but not yet calculated.
!!    * 2 if allocated and calculated
!!
!! SOURCE

 type,public :: melements_flags_type

  integer :: has_hbare=0
  integer :: has_lexexch=0
  integer :: has_sxcore=0
  integer :: has_vhartree=0
  integer :: has_vu=0
  integer :: has_vxc=0
  integer :: has_vxcval=0
  integer :: has_vxcval_hybrid=0
  integer :: only_diago=0
  ! 1 if only diagonal elements are calculated

 end type melements_flags_type
!!***

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

!!****t* m_melemts/melements_type
!! NAME
!!
!! FUNCTION
!!  Structure defining a database to store the matrix elements of operators
!!  needed for GW calculations.
!!
!! SOURCE

 type,public :: melements_type

  integer :: nkibz
  ! Number of k-points.

  integer :: nsppol 
  ! Number of independent spin-polarizations.

  integer :: nspinor 
  ! 1 for collinear, 2 for noncollinear.

  integer :: nspden
  ! Number of independent spin-density components.

  integer :: bmin,bmax
  ! min and Max band index over k-points and spin.
  ! Used to dimension the arrays below.
  
  integer, allocatable :: bands_idx(:,:,:)     
  ! bands_idx(2,nkibz,nsppol)
  ! min and Max band index for each k-point and spin.

  integer, allocatable :: iscalc(:,:)    
  ! stat_k(nkibz,nsppol)
  ! 1 if this k-point and spin has been calculated, 0 otherwise.

  real(dp), allocatable :: kibz(:,:)    
  ! kibz(3,nkibz)
  ! The list of k-points in reduced coordinates.

  complex(dpc), allocatable :: hbare(:,:,:,:)     
  ! hbare(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the bare Hamiltonian.

  complex(dpc), allocatable :: sxcore(:,:,:,:)   
  ! sxcore(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the Fock operator generated by core electrons.

  complex(dpc), allocatable :: vhartree(:,:,:,:)    
  ! vhartree(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the Hartree potential.

  complex(dpc), allocatable :: vlexx(:,:,:,:)    
  ! vlexx(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the local exact exchange potential.

  complex(dpc), allocatable :: vu(:,:,:,:)       
  ! vu(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the U Hamiltonian.

  complex(dpc), allocatable :: vxc(:,:,:,:)    
  ! vxc(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of XC potential, including core if present.

  complex(dpc), allocatable :: vxcval(:,:,:,:)    
  ! vxcval(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the XC potential, valence-only contribution.

  complex(dpc), allocatable :: vxcval_hybrid(:,:,:,:)    
  ! vxcval_hybrid(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the XC potential for hybrid calculations, valence-only contribution.

  type(melements_flags_type) :: flags

 end type melements_type
!!***

 integer,parameter,private :: NNAMES=8
 integer,parameter,private :: NAMELEN=13

! List of matrix element names.
! Important note: Remember to add a new entry here when you add a new matrix element.
! then change the routine my_select_melements accordingly.
 character(len=NAMELEN),parameter,private :: &
& ANAMES(NNAMES) = (/"vxc          ",&
&                    "vxcval       ",&
&                    "vxcval_hybrid",&
&                    "sxcore       ",&
&                    "vu           ",&
&                    "vlexx        ",&
&                    "vhartree     ",&
&                    "hbare        "/)

CONTAINS  !========================================================================================

!!****f* m_melemts/reset_mflags
!! NAME
!! reset_mflags
!!
!! FUNCTION
!!  Set all flags in melements_flags_type to 0.
!!
!! INPUTS
!!
!! PARENTS
!!      m_melemts,sigma
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine reset_mflags(Mflags)

 use defs_basis

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

 implicit none

!Arguments ------------------------------------
 type(melements_flags_type),intent(inout) :: Mflags

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

 ! @melements_flags_type
 Mflags%has_hbare           = 0
 Mflags%has_sxcore          = 0
 Mflags%has_vhartree        = 0 
 Mflags%has_vu              = 0
 Mflags%has_vxc             = 0
 Mflags%has_vxcval          = 0
 Mflags%has_vxcval_hybrid   = 0
 Mflags%has_lexexch         = 0
 Mflags%only_diago          = 0

end subroutine reset_mflags
!!***

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

!!****f* m_melemts/copy_mflags
!! NAME
!! copy_mflags
!!
!! FUNCTION
!!  Copy an object storing the flags.
!!
!! INPUTS
!!  Mflags_in=The flags to be copied.
!! 
!! OUTPUT
!!  Mflags_out=The new set of flags.
!!
!! PARENTS
!!      m_melemts
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine copy_mflags(Mflags_in, Mflags_out)

 use defs_basis

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

 implicit none

!Arguments ------------------------------------
 type(melements_flags_type),intent(in)    :: Mflags_in
 type(melements_flags_type),intent(inout) :: Mflags_out

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

 call reset_mflags(Mflags_out)

 ! @melements_flags_type
 Mflags_out%has_hbare           = Mflags_in%has_hbare  
 Mflags_out%has_sxcore          = Mflags_in%has_sxcore
 Mflags_out%has_vhartree        = Mflags_in%has_vhartree
 Mflags_out%has_vu              = Mflags_in%has_vu
 Mflags_out%has_vxc             = Mflags_in%has_vxc    
 Mflags_out%has_vxcval          = Mflags_in%has_vxcval 
 Mflags_out%has_vxcval_hybrid   = Mflags_in%has_vxcval_hybrid
 Mflags_out%has_lexexch         = Mflags_in%has_lexexch
 Mflags_out%only_diago          = Mflags_in%only_diago

end subroutine copy_mflags
!!***

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

!!****f* m_melemts/destroy_melements
!! NAME
!! destroy_melements
!!
!! FUNCTION
!!  Free all dynamic memory of the database
!!
!! INPUTS
!!  Mels<melements_type>=The database to be freed
!!
!! OUTPUT
!!  See side effects
!!
!! PARENTS
!!      sigma
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine destroy_melements(Mels)

 use defs_basis

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

 implicit none

!Arguments ------------------------------------
 type(melements_type),intent(inout) :: Mels

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

 DBG_ENTER("COLL")

 ! @melements_type

!integer arrays
 if (allocated(Mels%bands_idx)) then
   ABI_FREE(Mels%bands_idx)
 end if
 if (allocated(Mels%iscalc))   then
   ABI_FREE(Mels%iscalc)
 end if

!real arrays
 if (allocated(Mels%kibz)) then
   ABI_FREE(Mels%kibz)
 end if

!complex arrays
 if (allocated(Mels%hbare))   then
   ABI_FREE(Mels%hbare)
 end if
 if (allocated(Mels%sxcore)) then
   ABI_FREE(Mels%sxcore)
 end if
 if (allocated(Mels%vhartree))   then
   ABI_FREE(Mels%vhartree)
 end if
 if (allocated(Mels%vlexx)) then
   ABI_FREE(Mels%vlexx)
 end if
 if (allocated(Mels%vu)) then
   ABI_FREE(Mels%vu)
 end if
 if (allocated(Mels%vxc)) then
   ABI_FREE(Mels%vxc)
 end if
 if (allocated(Mels%vxcval)) then
   ABI_FREE(Mels%vxcval)
 end if
 if (allocated(Mels%vxcval_hybrid)) then
   ABI_FREE(Mels%vxcval_hybrid)
 end if

 ! * Reset all has_* flags.
 call reset_mflags(Mels%flags)

 DBG_EXIT("COLL")

end subroutine destroy_melements
!!***

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

!!****f* m_melemts/my_select_melements
!! NAME
!!  my_select_melements
!!
!! FUNCTION
!!  Helper function returning a pointer to the array "aname" as well as the status of the array
!!
!! INPUTS
!!  Mels<melements_type>=The database.
!!  aname=String with the name of the array.
!!
!! OUTPUT
!!  flag_p=Pointer to the integer defining the status of the array, see melements_flags_type.
!!  arr_p=The pointer to the array.
!!
!! PARENTS
!!      m_melemts
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine my_select_melements(Mels,aname,flag_p,arr_p)

 use defs_basis

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

 implicit none

!Arguments ------------------------------------
!scalars
 integer,pointer :: flag_p
 character(len=*),intent(in) :: aname
 type(melements_type),target,intent(in) :: Mels
!arrays
 complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)

!Local variables-------------------------------
 character(len=500) :: msg

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

 SELECT CASE (tolower(aname))
 CASE ("hbare") 
   flag_p => Mels%flags%has_hbare 
   arr_p  => Mels%hbare
 CASE ("sxcore")
   flag_p => Mels%flags%has_sxcore
   arr_p  => Mels%sxcore
 CASE ("vhartree") 
   flag_p => Mels%flags%has_vhartree
   arr_p  => Mels%vhartree
 CASE ("vlexx") 
   flag_p => Mels%flags%has_lexexch
   arr_p  => Mels%vlexx
 CASE ("vu") 
   flag_p => Mels%flags%has_vu
   arr_p  => Mels%vu
 CASE ("vxc") 
   flag_p => Mels%flags%has_vxc
   arr_p  => Mels%vxc
 CASE ("vxcval") 
   flag_p => Mels%flags%has_vxcval
   arr_p  => Mels%vxcval
 CASE ("vxcval_hybrid") 
   flag_p => Mels%flags%has_vxcval_hybrid
   arr_p  => Mels%vxcval_hybrid
 CASE DEFAULT
   msg = "Wrong aname = "//TRIM(aname)
   MSG_ERROR(msg)
 END SELECT 

end subroutine my_select_melements
!!***

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

!!****f* m_melemts/init_melements
!! NAME
!! init_melements
!!
!! FUNCTION
!!  Initialize the database, allocate arrays according to 
!!  Mflags_in, zeroing the content of the allocated arrays.
!!
!! INPUTS
!!  nspinor=Number of spinor components
!!  nsppol=Number of independent spin polarizations.
!!  nspden=Number of spin-density components
!!  b1,b2=min and Max band index over spin and k-points.
!!
!! OUTPUT
!!  Mels=The initialized database with dimensions and allocated memory.
!!
!! PARENTS
!!      calc_vhxc_me
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine init_melements(Mels,Mflags_in,nsppol,nspden,nspinor,nkibz,kibz,bands_idx)

 use defs_basis

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

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: nspinor,nspden,nsppol,nkibz
 type(melements_flags_type),intent(in) :: Mflags_in
 type(melements_type),intent(out) :: Mels
!arrays
 integer,intent(in) :: bands_idx(2,nkibz,nsppol)
 real(dp),intent(in) :: kibz(3,nkibz)

!Local variables-------------------------------
 integer :: ikibz,isppol,bmin,bmax,b1,b2

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

 DBG_ENTER("COLL")
 
 !@melements_type

 ! * Copy flags.
 call copy_mflags(Mflags_in, Mels%flags)

 ! * Copy dimensions.
 Mels%nkibz   = nkibz
 Mels%nsppol  = nsppol
 Mels%nspinor = nspinor
 Mels%nspden  = nspden

 ABI_MALLOC(Mels%bands_idx,(2,nkibz,nsppol))
 Mels%bands_idx=bands_idx

 ABI_MALLOC(Mels%iscalc,(nkibz,nsppol))
 Mels%iscalc=0

 bmin = HUGE(1); bmax =-HUGE(1)
 do isppol=1,Mels%nsppol
   do ikibz=1,Mels%nkibz
     if (ANY(Mels%bands_idx(:,ikibz,isppol)/=0)) then
       b1 = Mels%bands_idx(1,ikibz,isppol)
       b2 = Mels%bands_idx(2,ikibz,isppol)
       Mels%iscalc(ikibz,isppol)=1
       bmin = MIN(bmin,b1)
       bmax = MAX(bmax,b2)
       ABI_CHECK(b2>=b1 .and. b1>0,"Wrong b1, b2")
     end if
   end do
 end do

 if (bmin==HUGE(1).or.bmax==-HUGE(1)) then
   MSG_BUG("Wrong bands_idx")
 end if

 Mels%bmin = bmin
 Mels%bmax = bmax

 b1 = Mels%bmin      
 b2 = Mels%bmax      

! real arrays
 ABI_MALLOC(Mels%kibz,(3,nkibz))
 Mels%kibz = kibz

! complex arrays
 if (Mels%flags%has_hbare == 1) then 
   ABI_CALLOC(Mels%hbare,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_sxcore == 1) then 
   ABI_CALLOC(Mels%sxcore,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vhartree == 1) then 
   ABI_CALLOC(Mels%vhartree,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_lexexch == 1) then 
   ABI_CALLOC(Mels%vlexx,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vu == 1) then 
   ABI_CALLOC(Mels%vu,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vxc == 1) then 
   ABI_CALLOC(Mels%vxc,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vxcval == 1) then 
   ABI_CALLOC(Mels%vxcval,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vxcval_hybrid == 1) then 
   ABI_CALLOC(Mels%vxcval_hybrid,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 DBG_EXIT("COLL")

end subroutine init_melements
!!***

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

!!****f* m_melemts/herm_melements
!! NAME
!! herm_melements
!!
!! FUNCTION
!!  Reconstruc the lower triangle of all calculated arrays.
!!  Assuming Hermitian operator. Works both for collinear and 
!!  non-collinear cases.
!!
!! INPUTS
!!  Mels=The database 
!!  [aname]=The name of the array to be symmetrized, by default
!!    all calculated arrays are filled.
!!
!! SIDE EFFECTS
!!  All arrays whose flag is 2, are filled assuming an Hermitian operator.
!!
!! PARENTS
!!      calc_vhxc_me
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine herm_melements(Mels,aname)

 use defs_basis

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

 implicit none

!Arguments ------------------------------------
!scalars
 character(len=*),optional,intent(in) :: aname
 type(melements_type),intent(inout) :: Mels

!Local variables-------------------------------
 integer :: is,ik,ib,jb,iab,iab_tr,iname
 integer,pointer :: flag_p
 character(len=NAMELEN) :: key
!arrays 
 integer,parameter :: trsp_idx(2:4)=(/2,4,3/)
 complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)

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

 ! === Symmetrize matrix elements ===                                                                        
 ! * In the collinear case, generate the lower triangle by just doing a complex conjugate.
 ! * In the noncollinear case do also a transposition since A_{12}^{ab} = A_{21}^{ba}^*
 !   2-->2, 3-->4, 4-->3
 !
 do iname=1,NNAMES
   key = ANAMES(iname)
   if (PRESENT(aname)) then
     if (key /= aname) CYCLE
   end if
   
   call my_select_melements(Mels,key,flag_p,arr_p)

   if (flag_p>0) then
     do ik=1,Mels%nkibz 
       do is=1,Mels%nsppol

         do jb=Mels%bmin,Mels%bmax
           do ib=Mels%bmin,jb ! Upper triangle 
   
             if (ib/=jb) then 
               arr_p(jb,ib,ik,is)=CONJG(arr_p(ib,jb,ik,is))
               if (Mels%nspinor==2) then 
                 do iab=2,4
                   iab_tr=trsp_idx(iab)
                   arr_p(jb,ib,ik,iab)=CONJG(arr_p(ib,jb,ik,iab_tr))
                 end do
               end if
             else ! For ib==jb force real-valued
               arr_p(jb,ib,ik,is)=half*(arr_p(jb,ib,ik,is)+CONJG(arr_p(jb,ib,ik,is)))
               if (Mels%nspinor==2) arr_p(jb,ib,ik,2)=half*(arr_p(ib,jb,ik,2)+CONJG(arr_p(ib,jb,ik,2)))
             end if
   
           end do !ib
         end do !jb
   
       end do !is
     end do !ik
   end if

 end do !inames

end subroutine herm_melements
!!***

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

!!****f* m_melemts/xsum_melements
!! NAME
!! xsum_melements
!!
!! FUNCTION
!!  Perform a collective SUM within the MPI communicator comm 
!!  of the matrix elements stored in the database.
!!
!! INPUTS
!!  Mels=The database 
!!  [aname]=The name of a particular array to be summed, by default
!!    all allocated arrays are considered.
!!
!! SIDE EFFECTS
!!  All arrays whose flag==1 are summed within the MPI communicator comm.
!!  In output the corresponding flas is set to 2.
!!
!! PARENTS
!!      calc_vhxc_me
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine xsum_melements(Mels,comm,aname)

 use defs_basis

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

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: comm
 character(len=*),optional,intent(in) :: aname
 type(melements_type),intent(inout) :: Mels

!Local variables-------------------------------
 integer :: iname,ierr
 integer,pointer :: flag_p
 character(len=NAMELEN) :: key
 character(len=500) :: msg      
!arrays 
 complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)

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

 do iname=1,NNAMES
   key = ANAMES(iname)
   if (PRESENT(aname)) then
    if (key /= aname) CYCLE
   end if
   
   call my_select_melements(Mels,key,flag_p,arr_p)

   if (flag_p==1) then
     call xmpi_sum(arr_p,comm,ierr)
     if (ierr/=0) then
       write(msg,'(a,i4,2a)')" xmpi_sum reported ierr= ",ierr," for key= ",TRIM(key)
       MSG_ERROR(msg)
     end if
     flag_p=2 ! Tag this array as calculated
   end if
 end do

 call xmpi_barrier(comm)

end subroutine xsum_melements
!!***

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

!!****f* m_melemts/print_melements
!! NAME
!! print_melements
!!
!! FUNCTION
!!  Printout of the content of all calculated array.
!!  Optionally, it is possible to print the content of a 
!!  single entry of the database.
!!
!! INPUTS
!!  Mels=The database 
!!  [unit]=the unit number for output, defaults to std_out
!!  [prtvol]=verbosity level, defaults to 0
!!  [mode_paral]=either "COLL" or "PERS", default to "COLL"
!!  [header]=title for info
!!
!! OUTPUT
!!  Only writing
!!
!! PARENTS
!!      sigma
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine print_melements(Mels,names_list,header,unit,prtvol,mode_paral)

 use defs_basis

!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 'print_melements'
 use interfaces_14_hidewrite
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 type(melements_type),intent(in) :: Mels
 integer,optional,intent(in) :: prtvol,unit
 character(len=*),optional,intent(in) :: names_list(:)
 character(len=*),optional,intent(in) :: header
 character(len=4),optional,intent(in) :: mode_paral

!Local variables-------------------------------
 integer :: my_unt,my_prtvol,max_r,max_c,ii
 integer :: isppol,ikibz,iab,ib,b1,b2,my_nkeys,ikey
 integer,pointer :: flag_p
 character(len=4) :: my_mode
 character(len=NAMELEN) :: key
 character(len=500) :: msg,str,fmt
!arrays 
 integer,allocatable :: tab(:)
 character(len=NAMELEN),allocatable :: my_keys(:)
 complex(dpc),allocatable :: mat(:,:)

 type rarr_dpc4
   complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)
 end type rarr_dpc4
 type(rarr_dpc4),allocatable :: data_p(:)

! *************************************************************************
 
 !@melements_type
 my_unt   =std_out; if (PRESENT(unit      )) my_unt   =unit
 my_prtvol=0      ; if (PRESENT(prtvol    )) my_prtvol=prtvol
 my_mode  ='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral

 ABI_CHECK(Mels%nspinor==1,"nspinor=2 not coded")

 if (PRESENT(names_list)) then 
   my_nkeys=SIZE(names_list)
   ABI_MALLOC(my_keys,(my_nkeys))
   my_keys = names_list
 else 
   my_nkeys=NNAMES
   ABI_MALLOC(my_keys,(NNAMES))
   my_keys = ANAMES
 end if

 ABI_DT_MALLOC(data_p,(my_nkeys))
 ABI_MALLOC(tab,(my_nkeys))
 tab=0

 my_nkeys=0; str = "  ib"; ii=4
 do ikey=1,SIZE(my_keys)
   key = my_keys(ikey) 
   call my_select_melements(Mels,key,flag_p,data_p(ikey)%arr_p)
   if (flag_p==2) then 
     my_nkeys=my_nkeys+1
     tab(my_nkeys)=ikey
     str(ii+1:)=" "//TRIM(tolower(key))
     ii = ii+MAX(1+LEN_TRIM(key),10)
     ABI_CHECK(ii<490,"I'm gonna SIGFAULT!")
   end if
 end do
 
 write(msg,'(2a)')ch10,' === Matrix Elements stored in Mels% [eV] === '
 if (PRESENT(header)) write(msg,'(4a)')ch10,' === '//TRIM(ADJUSTL(header))//' [eV] === '
 call wrtout(my_unt,msg,my_mode)

 if (my_nkeys==0) GOTO 10
 write(fmt,'(a,i4,a)')'(1x,i3,',my_nkeys,'(1x,f9.5))' ! width of 10 chars

 do isppol=1,Mels%nsppol
   do ikibz=1,Mels%nkibz

    if (Mels%iscalc(ikibz,isppol)/=1) CYCLE
                                                                                  
    write(msg,'(a,3es16.8,a,i2,a)')" kpt= (",Mels%kibz(:,ikibz),") spin=",isppol,":"
    call wrtout(my_unt,msg,my_mode)

    b1 = Mels%bands_idx(1,ikibz,isppol)
    b2 = Mels%bands_idx(2,ikibz,isppol)

    if (Mels%flags%only_diago==1.or.my_prtvol==0) then ! Print only the diagonal.
      write(msg,'(a)')str
      call wrtout(my_unt,msg,my_mode)
      do ib=b1,b2
        do iab=1,Mels%nspinor**2
          !write(msg,'(i3,(f8.3))')ib,(REAL (arr_p(ib,ib,ikibz,iab))*Ha_eV, iab=1,Mels%nspinor**2)
          write(msg,fmt)ib,(REAL (data_p(tab(ikey))%arr_p(ib,ib,ikibz,iab))*Ha_eV, ikey=1,my_nkeys)
          call wrtout(my_unt,msg,my_mode)
        end do
      end do

    else ! Print full matrix.
      max_r = b2-b1+1
      max_c = MIN(b2-b1+1, 9)
      ABI_MALLOC(mat,(b1:b2,b1:b2))
      do ikey=1,my_nkeys
        write(msg,'(3a)')" **** Off-diagonal elements of ",TRIM(my_keys(tab(ikey)))," **** "
        call wrtout(my_unt,msg,my_mode)
        do iab=1,Mels%nspinor**2
          mat = data_p(tab(ikey))%arr_p(b1:b2,b1:b2,ikibz,iab)*Ha_eV
          call print_arr(mat,max_r,max_c,my_unt,my_mode)
        end do
        write(msg,'(a)')ch10
        call wrtout(my_unt,msg,my_mode)
      end do
      ABI_FREE(mat)
    end if

   end do !ikibz
 end do ! isppol

10 continue
 ABI_FREE(my_keys)
 ABI_DT_FREE(data_p)
 ABI_FREE(tab)

end subroutine print_melements
!!***

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

!!****f* m_melemts/zero_melements    
!! NAME
!! zero_melements
!!
!! FUNCTION
!!
!! INPUTS
!!   irrep_tab=Array used to select the entries that have to be set to zero. 
!!     irrep_tab(ib,ik,is)=gives the index of the irreducible representation associated to state (ib,ik,is).
!!  [aname]=The name of the array to be symmetrized, by default
!!    all calculated arrays are filled.
!!
!! SIDE EFFECTS
!!  Mels= All arrays elements connecting states belonging to different irreps are set to zero.
!!
!! PARENTS
!!      sigma
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine zero_melements(Mels,irrep_tab,aname)

 use defs_basis

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

 implicit none

!Arguments ------------------------------------
!scalars
 character(len=*),optional,intent(in) :: aname
 type(melements_type),intent(inout) :: Mels
!arrays
 integer,intent(in) :: irrep_tab(Mels%bmin:Mels%bmax,Mels%nkibz,Mels%nsppol)

!Local variables-------------------------------
 integer :: is,ik,ib,jb,iname,irrep_j,irrep_i
 integer,pointer :: flag_p
 character(len=NAMELEN) :: key
!arrays 
 complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)

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

 do iname=1,NNAMES
   key = ANAMES(iname)
   if (PRESENT(aname)) then
     if (key /= aname) CYCLE
   end if
   
   call my_select_melements(Mels,key,flag_p,arr_p)

   if (flag_p>0) then
     do is=1,Mels%nsppol
       do ik=1,Mels%nkibz 

         do jb=Mels%bmin,Mels%bmax
           irrep_j = irrep_tab(jb,ik,is)
           do ib=Mels%bmin,Mels%bmax
             irrep_i = irrep_tab(ib,ik,is)
             !
             ! Set this matrix element to zero if the irreps are known and they differ. 
             if (irrep_i/=irrep_j .and. ALL((/irrep_i,irrep_j/) /=0) ) then 
               !write(std_out,*)"setting to zero ",ib,jb,ik,is
               if (Mels%nspinor==2) then 
                 arr_p(ib,jb,ik,is)=czero
               else
                 arr_p(ib,jb,ik,:)=czero
               end if
             end if

           end do !ib
         end do !jb
   
       end do !is
     end do !ik
   end if

 end do !inames

end subroutine zero_melements
!!***

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

END MODULE m_melemts
!!***
