!{\src2tex{textfont=tt}}
!!****f* ABINIT/wfsinp
!! NAME
!! wfsinp
!!
!! FUNCTION
!! Do initialization of wavefunction files.
!! Also call other relevant routines for this initialisation.
!! Detailed description :
!!  - Initialize unit wff1 for input of wf data
!!
!! formeig option (format of the eigenvalues and occupations) :
!!   0 => ground-state format (initialisation of
!!        eigenvectors with random numbers, vector of eigenvalues,
!!        occupations are present)
!!   1 => respfn format (initialisation of
!!        eigenvectors with 0 s, hermitian matrix of eigenvalues)
!!
!! COPYRIGHT
!! Copyright (C) 1998-2014 ABINIT group (DCA, XG, GMR, AR)
!! 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
!!  ecut0=kinetic energy cutoffs for basis sphere 0 (hartree) (if squeeze=1)
!!  ecut=kinetic energy cutoffs beyond which the coefficients of cg vanish (Ha)
!!   (needed only if squeeze=1)
!!  ecut_eff=effective kinetic energy planewave cutoff (hartree), needed
!!    to generate the sphere of plane wave
!!  exchn2n3d=if 1, n2 and n3 are exchanged
!!  formeig=explained above
!!  gmet(3,3), gmet0(3,3)=reciprocal space metrics (bohr^-2)
!!  headform0=header format (might be needed to read the block of wfs)
!!  indkk(nkpt*sppoldbl,6)=describe k point number of kptns0 that allows to
!!   generate wavefunctions closest to given kpt
!!   indkk(:,1)=k point number of kptns0
!!   indkk(:,2)=symmetry operation to be applied to kpt0, to give kpt0a
!!    (if 0, means no symmetry operation, equivalent to identity )
!!   indkk(:,3:5)=shift in reciprocal space to be given to kpt0a,
!!    to give kpt0b, that is the closest to kpt.
!!   indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
!!  indkk0(nkpt0,nkassoc)=list of k points that will be generated by k point number ikpt0
!!  istwfk(nkpt)=input parameter that describes the storage of wfs
!!  istwfk0(nkpt0)=input parameter that describes the storage of wfs in set0
!!  kptns(3,nkpt),kptns0(3,nkpt0)=k point sets (reduced coordinates)
!!  localrdwf=(for parallel case) if 1, the wff1%unwff file is local to each machine
!!  mband=maximum number of bands
!!  mban_dp_rd=maximum number of bands to be read
!!  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
!!  mkmem=Number of k-points treated by this node.
!!  mpi_enreg=informations about MPI parallelization
!!  mpi_enreg0=informations about MPI parallelization in set0
!!  mpw=maximum number of planewaves as dimensioned in calling routine
!!  mpw0=maximum number of planewaves on disk file
!!  nban_dp_rd(nkpt0*nsppol0)=number of bands to be read at each k point
!!  nband(nkpt*nsppol)=number of bands at each k point
!!  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft
!!  nkassoc=dimension of indkk0 array
!!  nkpt=number of k points expected
!!  nkpt0=number of k points on disk
!!  npwarr(nkpt)=array holding npw for each k point.
!!  npwarr0(nkpt0)=array holding npw for each k point, disk format.
!!  nspinor=number of spinorial components of the wavefunctions
!!  nspinor0=number of spinorial components of the wavefunctions on disk
!!  nsppol=1 for unpolarized, 2 for spin-polarized
!!  nsppol0=1 for unpolarized, 2 for spin-polarized, when on disk
!!  nsym=number of symmetry elements in space group
!!  optorth= 1 if the WFS have to be orthogonalized; 0 otherwise
!!  prtvol=control print volume and debugging
!!  restart= if 2, want conversion between wavefunctions
!!           if 1, direct restart is allowed (see hdr_check.f)
!!  rprimd(3,3)=dimensional primitive translations for real space (bohr)
!!  sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn
!!            if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries
!!  squeeze=1 if cg_disk is to be used, even with mkmem/=0
!!  symrel(3,3,nsym)=symmetry operations in real space in terms
!!   of primitive translations
!!  tnons(3,nsym)=nonsymmorphic translations for symmetry operations
!!  wff1, structure informations for input and output files
!!
!! OUTPUT
!!  if ground state format (formeig=0):
!!    eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
!!  if respfn format (formeig=1):
!!    eigen(2*mband*mband*nkpt*nsppol)=
!!         matrix of eigenvalues (input or init to large number), (Ha)
!! Conditional output:
!!    cg_disk(2,mpw*nspinor*mband*mkmem*nsppol)=complex wf array
!!      be careful : an array of size cg(2,npw*nspinor), as used
!!      in the response function code, is not enough !
!!
!! SIDE EFFECTS
!!  if ground state format (formeig=0):
!!    occ(mband*nkpt*nsppol)=occupations (from disk or left at their initial value)
!!    NOT OUTPUT NOW !
!!
!! NOTES
!! occ will not be modified nor output, in the present status of this routine.
!!
!! WARNINGS
!! For parallelism : no distinction yet between nban_dp_rd and nband
!!
!! TODO
!! THE DESCRIPTION IS TO BE COMPLETELY REVISED, AS THIS ONE COMES FROM inwffil.f
!!
!! PARENTS
!!      inwffil
!!
!! CHILDREN
!!      initwf,mpi_bcast,mpi_recv,mpi_send,pareigocc,timab,wfconv,wffreadskipk
!!      wrtout
!!
!! SOURCE

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

#include "abi_common.h"


subroutine wfsinp(cg,cg_disk,ecut,ecut0,ecut_eff,eigen,exchn2n3d,&
&                  formeig,gmet,gmet0,headform0,indkk,indkk0,istwfk,&
&                  istwfk0,kptns,kptns0,localrdwf,mband,mban_dp_rd,&
&                  mcg,mcg_disk,mkmem,mpi_enreg,mpi_enreg0,mpw,mpw0,nband,nban_dp_rd,&
&                  ngfft,nkassoc,nkpt,nkpt0,npwarr,npwarr0,nspinor,&
&                  nspinor0,nsppol,nsppol0,nsym,occ,optorth,prtvol,restart,rprimd,&
&                  sppoldbl,squeeze,symrel,tnons,wff1)

 use defs_basis
 use defs_abitypes
 use m_errors
 use m_profiling
 use m_wffile
 use m_xmpi

#if defined HAVE_MPI2
 use mpi
#endif

!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 'wfsinp'
 use interfaces_14_hidewrite
 use interfaces_18_timing
 use interfaces_32_util
 use interfaces_62_iowfdenpot
 use interfaces_62_occeig
 use interfaces_66_wfs
!End of the abilint section

 implicit none

#if defined HAVE_MPI1
 include 'mpif.h'
#endif

!Arguments ------------------------------------
 integer, intent(in) :: exchn2n3d,formeig,headform0,localrdwf,mband,mban_dp_rd,mcg,mcg_disk,mkmem
 integer, intent(in) :: mpw,mpw0,nkassoc,nkpt,nkpt0,nspinor,nspinor0,nsppol,nsppol0,nsym
 integer, intent(in) :: optorth,prtvol,restart,sppoldbl,squeeze
 real(dp), intent(in) :: ecut,ecut0,ecut_eff
 type(MPI_type), intent(inout) :: mpi_enreg,mpi_enreg0
 type(wffile_type), intent(inout) :: wff1
 integer, intent(in) :: indkk(nkpt*sppoldbl,6),indkk0(nkpt0,nkassoc),istwfk(nkpt)
 integer, intent(in) :: istwfk0(nkpt0),nband(nkpt*nsppol),nban_dp_rd(nkpt0*nsppol0)
 integer, intent(in) :: ngfft(18),npwarr(nkpt),npwarr0(nkpt0),symrel(3,3,nsym)
 real(dp), intent(in) :: gmet(3,3),gmet0(3,3),kptns(3,nkpt),kptns0(3,nkpt0),rprimd(3,3)
 real(dp), intent(in) :: tnons(3,nsym)
 real(dp), intent(out) :: eigen((2*mband)**formeig*mband*nkpt*nsppol)
 real(dp), intent(inout) :: cg(2,mcg),cg_disk(2,mcg_disk) !vz_i pw_ortho
 real(dp), intent(inout) :: occ(mband*nkpt*nsppol)

!Local variables-------------------------------
 integer :: band_index,band_index_trial,ceksp,debug,dim_eig_k,iband,icg
 integer :: icg_disk,icg_trial,idum,ierr,ii,ikassoc,ikassoc_trial,ikpt,ikpt0
 integer :: ikpt10,ikpt_trial,ikptsp,ikptsp_old,inplace,isp,isp_max,isppol,isppol0
! integer :: ipw ! commented out below
 integer :: isppol_trial,me,mgfft,my_nspinor,my_nspinor0
 integer :: nban_dp_k,nban_dp_rdk,nband_k,nband_rdk,nband_trial,nbd,nbd_max
 integer :: ncopy,nkpt_eff,nproc_max,npw0_k,npw_k,npw_ktrial
 integer :: read_cg,read_cg_disk,sender,spaceComm
 character(len=500) :: message
 integer,allocatable :: band_index_k(:,:),icg_k(:,:),kg0_k(:,:),kg_k(:,:)
 real(dp) :: tsec(2)
 real(dp),allocatable :: eig0_k(:),eig_k(:),occ0_k(:),occ_k(:)
#if defined HAVE_MPI
 integer :: iproc,my_ikpt
 integer :: tag,test_cycle
 integer :: statux(MPI_STATUS_SIZE)
 integer,allocatable :: ikassoc_me(:),ikpt_me(:),isppol_me(:),nband_k_me(:)
#endif
 integer :: nkpt_max=50

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

!DEBUG
!write(std_out,*)' wfsinp : enter'
!write(std_out,*)' wfsinp : nband=',nband(:)
!write(std_out,*)' wfsinp : nban_dp_rd=',nban_dp_rd(:)
!write(std_out,*)' wfsinp : localrdwf=',localrdwf
!write(std_out,*)' wfsinp : paralbd,formeig=',mpi_enreg%paralbd,formeig
!write(std_out,*)' wfsinp : indkk0(:,1)=',indkk0(:,1)
!ENDDEBUG

 call timab(720,1,tsec)
 call timab(721,3,tsec)

 nkpt_max=50; if(xmpi_paral==1)nkpt_max=-1
 nbd_max=size(mpi_enreg%proc_distrb,2)
 isp_max=size(mpi_enreg%proc_distrb,3)

!Init mpi_comm
 spaceComm=mpi_enreg%comm_cell
 nproc_max=xcomm_size(spaceComm)
 me=mpi_enreg%me_kpt
 sender = 0

#if defined HAVE_MPI
 if(localrdwf==0)then
   ABI_ALLOCATE(ikpt_me,(nproc_max))
   ABI_ALLOCATE(nband_k_me,(nproc_max))
   ABI_ALLOCATE(ikassoc_me,(nproc_max))
   ABI_ALLOCATE(isppol_me,(nproc_max))
 end if
#endif

!Check the validity of formeig
 if(formeig/=0.and.formeig/=1)then
   write(message, '(a,i0,a)' )' formeig=',formeig,' , but the only allowed values are 0 or 1.'
   MSG_BUG(message)
 end if

 my_nspinor =max(1,nspinor /mpi_enreg%nproc_spinor)
 my_nspinor0=max(1,nspinor0/mpi_enreg%nproc_spinor)

 nkpt_eff=max(nkpt0,nkpt)
 if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max

 ABI_ALLOCATE(icg_k,(nkpt,nsppol))
 ABI_ALLOCATE(band_index_k,(nkpt,nsppol))

!write(std_out,*)' wfsinp : me,isppol,ikpt,icg_k,band_index_k'

!Compute the locations of the blocks in cg, eig and occ
 icg=0
 band_index=0
 ikpt10=0

 do isppol=1,nsppol
   do ikpt=1,nkpt

     nband_k=nband(ikpt+(isppol-1)*nkpt)
     band_index_k(ikpt,isppol)=band_index

#if defined HAVE_MPI
     test_cycle=0;nbd=min(nband_k,nbd_max);isp=min(isppol,isp_max)
     if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nbd,isp,me))test_cycle=1
     if(test_cycle==1)then
       band_index=band_index+nband_k*(2*nband_k)**formeig
!      In the case this k point does not belong to me, cycle
       cycle
     end if
#endif

     npw_k=npwarr(ikpt)
     icg_k(ikpt,isppol)=icg
     icg=icg+npw_k*my_nspinor*nband_k

     band_index=band_index+nband_k*(2*nband_k)**formeig
!    write(std_out,'(5i8)' )me,isppol,ikpt,icg_k(ikpt,isppol),band_index_k(ikpt,isppol)
   end do ! End k point loop
 end do  ! End spin loop

 band_index=0
 ikptsp_old=0

!DEBUG
!write(std_out,*)' wfsinp: before loop'
!write(std_out,*)' nsppol0,nsppol,nkpt0',nsppol0,nsppol,nkpt0
!write(std_out,*)' mpw,mgfft,mpw,mpw0',mpw,mgfft,mpw,mpw0
!ENDDEBUG

 mgfft=maxval(ngfft(1:3))
 if(squeeze==1)then
   ABI_ALLOCATE(kg_k,(3,mpw))
   ABI_ALLOCATE(kg0_k,(3,mpw0))
 end if

 eigen(:)=0.0_dp
!occ(:)=0.0_dp

 call timab(721,2,tsec)

!Loop over spins
!For the time being, do not allow nsppol=2 to nspinor=2 conversion
!MT 20110707: this can be done by a fake call to the routine: see inwffil
 do isppol0=1,min(nsppol0,nsppol)

!  Loop on k points  : get the cg then eventually write on unwfnow
   do ikpt0=1,nkpt0

     call timab(722,1,tsec)

     nban_dp_rdk=nban_dp_rd(ikpt0+(isppol0-1)*nkpt0)

!    DEBUG
!    write(std_out,*)' wfsinp: ikpt0,isppol0,nkpt0=',ikpt0,isppol0,nkpt0
!    write(std_out,*)' nban_dp_rdk=',nban_dp_rdk
!    ENDDEBUG

     npw0_k=npwarr0(ikpt0)
     if(ikpt0<=nkpt_eff)then
       write(message,'(a,a,2i4)')ch10,' wfsinp: inside loop, init ikpt0,isppol0=',ikpt0,isppol0
       call wrtout(std_out,message,'PERS')
     end if

!    Must know whether this k point is needed, and in which
!    block (ikpt, isppol), the wavefunction is to be placed.
!    Select the one for which the number of bands is the biggest.
     ikpt=0
     isppol=0
     ikassoc=0
     nband_k=0
#if defined HAVE_MPI
     if(localrdwf==0)then
       nband_k_me(:)=0
       ikpt_me(:)=0
       isppol_me(:)=0
       ikassoc_me(:)=0
       nband_k_me(:)=0
     end if
#endif

     do isppol_trial=1,nsppol

       if(nsppol==2 .and. nsppol0==2 .and. isppol0/=isppol_trial)cycle

       do ikassoc_trial=1,nkassoc

         ikpt_trial=indkk0(ikpt0,ikassoc_trial)
         if(sppoldbl==2)then
           if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle
           if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle
           if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt
         end if

#if defined HAVE_MPI
         if(localrdwf==1)then
           if(ikpt_trial/=0)then
             nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
             nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
             if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me))ikpt_trial=0
           end if
         end if
#endif

         if(ikpt_trial/=0)then
           nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
           if(nband_k<nband_trial)then
             nband_k=nband_trial ; ikpt=ikpt_trial ; isppol=isppol_trial
             ikassoc=ikassoc_trial
           end if

#if defined HAVE_MPI
           if(localrdwf==0)then
             do iproc=1,nproc_max
               my_ikpt=1
               nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
               nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
               if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,(iproc-1))) my_ikpt=0
               if(my_ikpt/=0)then
                 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
                 if(nband_k_me(iproc)<nband_trial)then
                   nband_k_me(iproc)=nband_trial
                   ikpt_me(iproc)=ikpt_trial
                   isppol_me(iproc)=isppol_trial
                   ikassoc_me(iproc)=ikassoc_trial
                 end if
               end if
             end do
           end if
#endif

         end if

       end do ! ikassoc_trial
     end do ! isppol_trial

!    DEBUG
!    write(std_out,*)' wfsinp : me,select isppol,ikpt=',me,isppol,ikpt
#if defined HAVE_MPI
!    write(std_out,*)' wfsinp : me,ikpt_me(:)=',me,ikpt_me(:)
#endif
!    write(std_out,*)' wfsinp : me,isppol_me(:)=',me,isppol_me(:)
!    stop
!    ENDDEBUG

     call timab(722,2,tsec)

!    If the wavefunction block to be read is interesting ...
     if (ikpt/=0)then
       
       call timab(723,3,tsec)
       sender = me
       nband_k=nband(ikpt+(isppol-1)*nkpt)
       npw_k=npwarr(ikpt)

#if defined HAVE_MPI
       if (localrdwf==1.or.(localrdwf==0.and.me==0)) then

         if(ikpt<=nkpt_eff)then
           write(message,'(a,i6,a,i8,a,i4,a,i4)') &
&           ' wfsinp: treating ',nband_k,' bands with npw=',npw_k,' for ikpt=',ikpt,' by node ',me
           call wrtout(std_out,message,'PERS')
         else if(ikpt==nkpt_eff+1)then
           call wrtout(std_out,' wfsinp: prtvol=0 or 1, do not print more k-points.','PERS')
         end if

       end if
#endif

       nband_rdk=nban_dp_rdk
       if(squeeze==1)nband_rdk=(nban_dp_rdk/nspinor0)*nspinor
       if(formeig==0)then
         ABI_ALLOCATE(eig_k,(nband_rdk))
         ABI_ALLOCATE(occ_k,(nband_rdk))
         ABI_ALLOCATE(eig0_k,(nban_dp_rdk))
         ABI_ALLOCATE(occ0_k,(nban_dp_rdk))
         dim_eig_k=nband_rdk
       else if(formeig==1)then
         ABI_ALLOCATE(eig_k,(2*nband_rdk*nband_rdk))
         ABI_ALLOCATE(eig0_k,(2*nban_dp_rdk*nban_dp_rdk))
         dim_eig_k=2*nband_rdk*nband_rdk
         ABI_ALLOCATE(occ0_k,(0))
         ABI_ALLOCATE(occ_k,(0))
       else
         ABI_ALLOCATE(occ0_k,(0))
         ABI_ALLOCATE(eig0_k,(0))
         ABI_ALLOCATE(occ_k,(0))
         ABI_ALLOCATE(eig_k,(0))
       end if
       eig_k(:)=0.0_dp
       eig0_k(:)=0.0_dp

!      Generate or read the cg for this k point
!      Either read into cg, or read into cg_disk
       read_cg=1 ; read_cg_disk=0
       if(squeeze==1)then
         read_cg=0 ; read_cg_disk=1
       end if
#if defined HAVE_MPI
       if(localrdwf==0)then
         read_cg=0
         read_cg_disk=0
!        XG20040106 The following condition is correct
         if(me==0)read_cg_disk=1
       end if
#endif

       icg=0
       if(read_cg==1)icg=icg_k(ikpt,isppol)

!      DEBUG
!      write(std_out,*)' wfsinp: before initwf',wff1%offwff
!      write(std_out,*)' wfsinp: me,read_cg,read_cg_disk=',me,read_cg,read_cg_disk
!      write(std_out,*)' wfsinp: nban_dp_rdk=',nban_dp_rdk
!      ENDDEBUG
       call timab(723,2,tsec)
       call timab(724,3,tsec)

       if(read_cg_disk==1)then
         call initwf (cg_disk,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,&
&         isppol0,mcg_disk,mpi_enreg0, &
&         nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1)
       end if

       if(read_cg==1)then
         call initwf (cg,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,&
&         isppol0,mcg,mpi_enreg0,&
&         nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1)
       end if

       call timab(724,2,tsec)
       call timab(725,3,tsec)

       nban_dp_k=min(nban_dp_rdk,(nband_k/nspinor)*nspinor0)
!      This band_index is defined BEFORE the eventual redefinition
!      of ikpt and isppol, needed  when localrdwf==0 in parallel
       band_index=band_index_k(ikpt,isppol)
!      DEBUG
!      write(std_out,*)' wfsinp: me,cg_disk(:,1)=',me,cg_disk(:,1)
!      ENDDEBUG

!      DEBUG
!      if(me==0 .and. ikpt0==1)then
!      write(std_out,*)' wfsinp : cg array, before trial, ikpt0=',ikpt0
!      do ipw=1,15
!      write(std_out,'(i4,2es20.10)' )ipw,cg(:,ipw)
!      end do
!      end if
!      ENDDEBUG


#if defined HAVE_MPI
       if(localrdwf==0)then
!        Warning: In that case , not yet // on nspinors
!        Transmit to each of the other processors, when needed
         if(nproc_max>=2)then
           do iproc=2,nproc_max
!            Only me=0 and me=iproc-1 are concerned by this
             if(me==0 .or. me==iproc-1)then

               ikpt=ikpt_me(iproc)
               isppol=isppol_me(iproc)

               if(ikpt/=0)then
!                In this case, processor iproc-1 needs the data
!                Generate a common tag
                 tag=256*(ikpt-1)+iproc+1
                 if(isppol==2)tag=-tag
                 nband_k=nband(ikpt+(isppol-1)*nkpt)
                 npw_k=npwarr(ikpt)
!                SEND
                 if(me==0)then
                   write(std_out,*)'SENDWFSINP ',me
                   call MPI_SEND(cg_disk,2*npw_k*my_nspinor*nband_k,&
&                   MPI_DOUBLE_PRECISION,iproc-1,tag,spaceComm,ierr)
                 end if
!                RECEIVE
                 if(me==iproc-1)then
                   call MPI_RECV(cg_disk,2*npw_k*my_nspinor*nband_k,&
&                   MPI_DOUBLE_PRECISION,0,tag,spaceComm,statux,ierr)
                   icg=icg_k(ikpt,isppol)
                   if(squeeze==0)then
                     cg(:,icg+1:icg+npw_k*my_nspinor*nband_k)=&
&                     cg_disk(:,1:npw_k*my_nspinor*nband_k)
                   end if
                   ikassoc=ikassoc_me(iproc)
                 end if
               end if

             end if
           end do ! iproc
         end if

!        DEBUG
!        write(std_out,*)' wfsinp: me, iproc loop finished',me
!        ENDDEBUG

!        Take care of me=0 needing the data
         if (me==0) then
           ikpt=ikpt_me(me+1)
           isppol=isppol_me(me+1)
           if(ikpt/=0 )then
             nband_k=nband(ikpt+(isppol-1)*nkpt)
             npw_k=npwarr(ikpt)
!            I am the master node, and I might need my own data
             icg=icg_k(ikpt,isppol)
             if(squeeze==0)then
!              Copy from cg_disk to cg
               cg(:,1+icg:npw_k*my_nspinor*nband_k+icg)= &
&               cg_disk(:,1:npw_k*my_nspinor*nband_k)
             end if
             ikassoc=ikassoc_me(me+1)
           end if
         end if
!        For the eigenvalues and occ, the transmission is much easier to write !
         call MPI_BCAST(eig0_k,nban_dp_rdk*(2*nban_dp_rdk)**formeig ,&
&         MPI_DOUBLE_PRECISION,0,spaceComm,ierr)
       end if
#endif

       if(formeig==0)then
!        The transfer from eig0_k to eig_k uses nban_dp_rdk, which contains
!        the maximal information.
         if(nspinor0==nspinor .or. squeeze==0)then
           eig_k(1:nban_dp_rdk)=eig0_k(1:nban_dp_rdk)
           occ_k(1:nban_dp_rdk)=occ0_k(1:nban_dp_rdk)
         else if(nspinor0==1 .and. nspinor==2)then
           do iband=1,nban_dp_rdk
             eig_k(2*iband  )=eig0_k(iband)
             eig_k(2*iband-1)=eig0_k(iband)
             occ_k(2*iband  )=occ0_k(iband)*0.5_dp
             occ_k(2*iband-1)=occ0_k(iband)*0.5_dp
           end do
         else if(nspinor0==2 .and. nspinor==1)then
           do iband=1,nban_dp_rdk
             eig_k(iband)=eig0_k(2*iband-1)
             occ_k(iband)=occ0_k(2*iband-1)*2.0_dp
           end do
         end if

!        DEBUG
!        write(std_out,*)' wfsinp: me,band_index,ikpt,isppol',me,band_index,ikpt,isppol
!        ENDDEBUG

!        The transfer to eigen uses nban_dp_k, that is bound by the number
!        of bands for this k point.
         ncopy=min(dim_eig_k,(nban_dp_k/nspinor0)*nspinor)
         eigen(1+band_index:ncopy+band_index)=eig_k(1:ncopy)
!        The transfer of occ should be done.

#if defined HAVE_MPI
         if(localrdwf==0 .and. ikpt/=0)then
!          Do not forget : ikpt,isppol were redefined ...
           band_index=band_index_k(ikpt,isppol)
           eigen(1+band_index:(nban_dp_k/nspinor0)*nspinor+band_index) = eig_k(1:(nban_dp_k/nspinor0)*nspinor)
!          The transfer of occ should be done.
         end if
#endif

       else if(formeig==1)then
         call wrtout(std_out,ABI_FUNC//': transfer of first-order eigs not yet coded!',"COLL")
       end if

!      DEBUG
!      write(std_out,*)' wfsinp : me,transferred eig_k',me
!      write(std_out,*)' me,mkmem,nsppol,nsppol0,isppol',me,mkmem,nsppol,nsppol0,isppol
!      write(std_out,*)' me,nkassoc,ikassoc',me,nkassoc,ikassoc
!      ENDDEBUG

       call timab(725,2,tsec)

!      Write to disk if appropriate
!      The coding has to be done ... here, only fragments ...
       call timab(727,3,tsec)

       do isppol_trial=1,nsppol
         if(nsppol==2 .and. nsppol0==2 .and. isppol_trial/=isppol)cycle
         do ikassoc_trial=1,nkassoc

!            DEBUG
!            write(std_out,*)' wfsinp: me, for ikassoc,isppol',&
!            &       me,ikassoc,isppol
!            write(std_out,*)' wfsinp: me, try ikassoc_trial,isppol_trial,nband_k',&
!            &       me,ikassoc_trial,isppol_trial,nband_k
!            ENDDEBUG

!            No conversion is to be done : it will be converted in newkpt
!            If squeeze==0, the block with the ikpt corresponding to ikassoc,
!            and with isppol, contains the wavefunction already
           if( ikassoc_trial/=ikassoc                   .or. &
&           (isppol_trial/=isppol .and. sppoldbl==1 ).or. &
&           squeeze==1                                       )then

             ikpt_trial=indkk0(ikpt0,ikassoc_trial)
             if(sppoldbl==2)then
               if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle
               if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle
               if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt
             end if


#if defined HAVE_MPI
             if(ikpt_trial/=0)then
               nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
               nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
               if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me)) ikpt_trial=0
             end if
#endif

             if(ikpt_trial/=0 .and. ikpt_trial<=nkpt_eff)then
               write(message,'(2a,2i5)')ch10,' wfsinp: transfer to ikpt_trial,isppol_trial=',ikpt_trial,isppol_trial
               call wrtout(std_out,message,'PERS')
             end if

             if(ikpt_trial/=0)then
               icg_trial=icg_k(ikpt_trial,isppol_trial)
               band_index_trial=band_index_k(ikpt_trial,isppol_trial)
               nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
               nban_dp_k=min(nban_dp_rdk,(nband_trial/nspinor)*nspinor0)

               if(squeeze==0)then
!                  GMR: modified to avoid compiler bug
!                  cg(:,1+icg_trial:npw0_k*my_nspinor0*nband_trial+icg_trial)=&
!                  &          cg(:,1+icg:npw0_k*my_nspinor0*nband_trial+icg)
                 do ii=1,npw0_k*my_nspinor0*nband_trial
                   cg(:,ii+icg_trial)=cg(:,ii+icg)
                 end do
!                  GMR
                 if(formeig==0)then
                   eigen(1+band_index_trial:nban_dp_k+band_index_trial)=eig_k(1:nban_dp_k)
!                    occ(1+band_index_trial:nban_dp_k+band_index_trial)=&
!                    &           occ_k(1:nban_dp_k)
                 end if
!                  RF transfer of eigenvalues still to be coded
               else if(squeeze==1)then
                 npw_ktrial=npwarr(ikpt_trial)
                 nband_k=(nban_dp_k/nspinor0)*nspinor
!                  Conversion to be done
                 ceksp=0 ; debug=0 ; icg_disk=0 ; idum=0 ; inplace=0
!                  Note that this routine also convert eig and occ
!                  even if the conversion had already been done


                 call wfconv(ceksp,cg_disk,cg,debug,ecut0,ecut,ecut_eff,&
&                 eig0_k,eig_k,exchn2n3d,formeig,gmet0,gmet,&
&                 icg_disk,icg_trial,ikpt0,ikpt10,ikpt_trial,indkk,&
&                 inplace,isppol_trial,istwfk0,istwfk,&
&                 kg0_k,kg_k,kptns0,kptns,nban_dp_rdk,nband_rdk,&
&                 mcg_disk,mcg,mpi_enreg0,mpi_enreg,mpw0,mpw,&
&                 nban_dp_rdk,nband_trial,ngfft,ngfft,nkpt0,nkpt,&
&                 npw0_k,npw_ktrial,nspinor0,nspinor,nsym,&
&                 occ0_k,occ_k,optorth,restart,rprimd,&
&                 sppoldbl,symrel,tnons)

!                  DEBUG
!                  write(std_out,*)' wfsinp: ikpt_trial=',ikpt_trial
!                  write(std_out,*)' nband_k,band_index_trial',nband_k,band_index_trial
!                  write(std_out,*)eig0_k(1:nban_dp_k)
!                  write(std_out,*)eig_k(1:nband_k)
!                  ENDDEBUG
                 eigen(1+band_index_trial:nband_k+band_index_trial)=eig_k(1:nband_k)
!                  occ(1+band_index_trial:nband_k+band_index_trial)=&
!                  &           occ_k(1:nband_k)
!                  RF transfer of eigenvalues still to be coded

!                  Endif squeeze==1
               end if

!                DEBUG
!                if(ikpt_trial==2)then
!                write(std_out,*)' wfsinp: iband,ipw,cg for ikpt_trial=2'
!                write(std_out,*)' nband_trial,npw_ktrial=',nband_trial,npw_ktrial
!                do iband=1,nband_trial
!                do ipw=1,npw_ktrial
!                write(std_out,'(2i5,2es16.6)' )&
!                &             iband,ipw,cg(:,ipw+(iband-1)*npw_ktrial+icg_trial)
!                end do
!                end do
!                end if
!                ENDDEBUG

!                End if ikpt_trial/=0
             end if

!              End if ikpt_trial already initialized
           end if

         end do ! ikassoc_trial
       end do ! isppol_trial

       call timab(727,2,tsec)

       ABI_DEALLOCATE(eig_k)
       ABI_DEALLOCATE(eig0_k)
       !if(formeig==0) then
         ABI_DEALLOCATE(occ_k)
         ABI_DEALLOCATE(occ0_k)
       !end if

     end if  ! End the condition of need of this k point
   end do ! End of the k loop
 end do !  End of spin loop

#if defined HAVE_MPI
 call timab(67,1,tsec)

!Still need to skip last k points to read history (MS)
!WARNING : not yet for formeig=1, but is it needed ?
 if(formeig==0)then
   do ikptsp=ikptsp_old+1,nkpt0*nsppol0
     isppol=1 ; if(ikptsp>nkpt0)isppol=2
     ikpt=ikptsp-nkpt0*(isppol-1)
     call WffReadSkipK(formeig,headform0,ikpt,isppol,mpi_enreg,wff1)
   end do
 end if

!Transmit eigenvalues. This routine works in both localrdwf=0 or 1 cases.
 call pareigocc(eigen,formeig,localrdwf,mpi_enreg,mband,nband,nkpt,nsppol,occ,1)


 if(localrdwf==0)then
   ABI_DEALLOCATE(ikpt_me)
   ABI_DEALLOCATE(nband_k_me)
   ABI_DEALLOCATE(ikassoc_me)
   ABI_DEALLOCATE(isppol_me)
 end if

 call timab(67,2,tsec)
#endif

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

 if(squeeze==1)then
   ABI_DEALLOCATE(kg_k)
   ABI_DEALLOCATE(kg0_k)
 end if


!DEBUG
!if(me==0)then
!write(std_out,*)' wfsinp : cg array='
!icg=0
!do isppol=1,nsppol
!do ikpt=1,1
!nband_k=nband(ikpt+(isppol-1)*nkpt)
!npw_k=npwarr(ikpt)
!do iband=1,nband_k
!write(std_out,*)' new band, icg=',icg
!do ipw=1,npw_k
!write(std_out,'(4i4,2es20.10)' )isppol,ikpt,iband,ipw,cg(:,icg+ipw)
!end do
!icg=icg+npw_k
!end do
!end do
!end do
!end if
!if(ireadwf==1)stop
!write(std_out,*)' wfsinp : eigen array='
!do ikpt=1,nkpt
!do iband=1,mband
!write(std_out,*)'ikpt,iband,eigen',ikpt,iband,eigen(iband+(ikpt-1)*mband)
!end do
!end do
!ENDDEBUG

 ABI_DEALLOCATE(icg_k)
 ABI_DEALLOCATE(band_index_k)

 call timab(720,2,tsec)

end subroutine wfsinp

!!***
