Logo Search packages:      
Sourcecode: abinit version File versions  Download package

wfsinp.F90

!{\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
!!  - Opens file on unit wffnow
!!    if the storage on disk is needed (mkmem==0)
!!    - Initializes wf data on wffnow, by calling the appropriate routine.
!!
!! 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)
!!
!! The wavefunctions after this initialisation are stored in unit wffnow%unwff
!!
!! COPYRIGHT
!! Copyright (C) 1998-2007 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
!!  exchn2n3=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
!!  mkmem=maximum number of k-points in core memory
!!  mpi_enreg=informations about MPI parallelization
!!  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
!!  symafm(nsym)=(anti)ferromagnetic part of symmetry operations
!!  symrel(3,3,nsym)=symmetry operations in real space in terms
!!   of primitive translations
!!  tnons(3,nsym)=nonsymmorphic translations for symmetry operations
!!  wff1,wffnow= 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:
!!    if mkmem==0:
!!    cg(2,mcg_disk)=partial complex wf array
!!                   allocated only when cg uses no space
!!    if mkmem/=0:
!!    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,leave_new,leave_test,mpi_bcast,mpi_recv,mpi_send,pareigocc,rwwf
!!      timab,wfconv,wffreadskipk,wrtout,xcomm_init,xmaster_init,xme_init
!!      xproc_init
!!
!! SOURCE

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

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

 use defs_basis
 use defs_datatypes

!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifdef HAVE_FORTRAN_INTERFACES
 use interfaces_01manage_mpi
 use interfaces_13io_mpi
 use interfaces_14iowfdenpot
 use interfaces_14wfs
 use interfaces_18seqpar, except_this_one => wfsinp
 use interfaces_lib01hidempi
#endif
!End of the abilint section

 implicit none

#if defined MPI
           include 'mpif.h'
#endif
!Arguments ------------------------------------
 integer, intent(in) :: exchn2n3,formeig,headform0,localrdwf,mband,mban_dp_rd,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
 type(wffile_type), intent(inout) :: wffnow, 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),symafm(nsym),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(out) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),cg_disk(2,mcg_disk)
 real(dp), intent(inout) :: occ(mband*nkpt*nsppol)

!Local variables-------------------------------
 integer :: band_index,band_index_trial,ceksp,debug,dim_eig_k,iband,iban_dp,icg
 integer :: icg_disk,icg_trial,idum,ierr,ii,ikassoc,ikassoc_trial,ikpt,ikpt0
 integer :: ikpt10,ikpt_trial,ikptsp,ikptsp_old,inplace,ipw,isppol,isppol0
 integer :: isppol_trial,isym,master,mcg,me,mgfft,nban_dp_k,nban_dp_rdk,nband_k
 integer :: nband_rdk,nband_trial,ncopy,nkpt_eff,nproc_max,npw0_k,npw_k
 integer :: npw_ktrial,read_cg,read_cg_disk,sender,spaceComm,tim_rwwf
 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(:)
!no_abirules
#if defined MPI
          !Variables introduced for MPI version
           integer :: case_transmit,iproc,my_ikpt,nbanddum,npwdum,nspinordum
           integer :: tag,test_cycle
           integer :: status(MPI_STATUS_SIZE)
           integer,allocatable :: ikassoc_me(:),ikpt_me(:),isppol_me(:),nband_k_me(:)
           real(dp),allocatable :: cg_send(:,:)
#endif
#if !defined MPI
           integer,parameter :: nkpt_max=50
#endif
#if defined MPI
           integer,parameter :: nkpt_max=-1
#endif

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

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

!Init me
 call xme_init(mpi_enreg,me)
!Init master
 call xmaster_init(mpi_enreg,master)
!Define nproc_max
 call xproc_init(mpi_enreg,nproc_max)
!Init mpi_comm
 call xcomm_init(mpi_enreg,spaceComm)
 sender = master

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

!Check the validity of formeig
 if(formeig/=0.and.formeig/=1)then
  write(message, '(a,a,a,a,i12,a)' ) ch10,&
&   ' wfsinp: BUG -',ch10,&
&   '  formeig=',formeig,' , but the only allowed values are 0 or 1.'
  call wrtout(06,message,'COLL')
  call leave_new('COLL')
 end if

 mcg=mpw*nspinor*mband*mkmem*nsppol

 if(mkmem==0)then
  write(6,*)' wfsinp : mkmem==0 not yet coded, stop'
  stop
 end if

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

 if(mkmem/=0)allocate(icg_k(nkpt,nsppol))
 allocate(band_index_k(nkpt,nsppol))

!DEBUG
!write(6,*)' wfsinp : me,isppol,ikpt,icg_k,band_index_k'
!stop
!ENDDEBUG

!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 MPI
          test_cycle=0
        if (mpi_enreg%parareel == 0) then
          if(minval(abs(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol)-me))&
&                /=0)test_cycle=1
        else
          if(mpi_enreg%proc_distrb_para(mpi_enreg%ipara,ikpt) /= me)&
&                test_cycle=1
        end if
          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)
   if(mkmem/=0)icg_k(ikpt,isppol)=icg
   icg=icg+npw_k*nspinor*nband_k
   band_index=band_index+nband_k*(2*nband_k)**formeig
!DEBUG
!   write(6,'(5i8)' )&
!&   me,isppol,ikpt,icg_k(ikpt,isppol),band_index_k(ikpt,isppol)
!ENDDEBUG

! End k point loop
  end do

!End spin loop
 end do

 band_index=0
 ikptsp_old=0

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

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

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

!Loop over spins
!For the time being, do not allow nsppol=2 to nspinor=2 conversion
 do isppol0=1,min(nsppol0,nsppol)

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

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

!DEBUG
!  write(6,*)' wfsinp: ikpt0,isppol0,nkpt0=',ikpt0,isppol0,nkpt0
!  write(6,*)' 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(6,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 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 MPI
           if(localrdwf==1)then
            if(ikpt_trial/=0)then
              nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
        if (mpi_enreg%parareel == 0) then
              if(minval(abs(mpi_enreg%proc_distrb(ikpt_trial,1:nband_trial,&
&                isppol_trial)-me))/=0) ikpt_trial=0
        else
              if(mpi_enreg%proc_distrb_para(mpi_enreg%ipara,ikpt_trial) /= me ) &
&                 ikpt_trial=0
        end if
            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 MPI
           if(localrdwf==0)then
            do iproc=1,nproc_max
             my_ikpt=1
             nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
        if (mpi_enreg%parareel == 0) then
             if(minval(abs(mpi_enreg%proc_distrb(ikpt_trial,1:nband_trial,&
&                isppol_trial)-(iproc-1)))/=0)  my_ikpt=0
        else
             if(mpi_enreg%proc_distrb_para(mpi_enreg%ipara,ikpt_trial) &
&                /= (iproc-1))  my_ikpt=0
        end if
             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(6,*)' wfsinp : me,select isppol,ikpt=',me,isppol,ikpt
#if defined MPI
!          write(6,*)' wfsinp : me,ikpt_me(:)=',me,ikpt_me(:)
#endif
!  write(6,*)' wfsinp : me,isppol_me(:)=',me,isppol_me(:)
!  stop
!ENDDEBUG

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

#if defined MPI
           if( localrdwf==1 .or. &
&             (localrdwf==0 .and.              &
&                  ( (me==0 .and. mpi_enreg%parareel == 0) .or. &
&                    (mpi_enreg%me_group_para==0 .and. mpi_enreg%parareel == 1)) &
&              )&
&            ) 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(06,message,'PERS')
            else if(ikpt==nkpt_eff+1)then
             write(message, '(a)' ) &
          &   ' wfsinp: prtvol=0 or 1, do not print more k-points.'
             call wrtout(06,message,'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
     allocate(eig_k(nband_rdk),occ_k(nband_rdk))
     allocate(eig0_k(nban_dp_rdk),occ0_k(nban_dp_rdk))
     dim_eig_k=nband_rdk
    else if(formeig==1)then
     allocate(eig_k(2*nband_rdk*nband_rdk))
     allocate(eig0_k(2*nban_dp_rdk*nban_dp_rdk))
     dim_eig_k=2*nband_rdk*nband_rdk
    end if
    eig_k(:)=0.0_dp
    eig0_k(:)=0.0_dp

!DEBUG
!   write(6,*)' me,ikpt,isppol,icg=',me,ikpt,isppol,icg
!   stop
!   if(ikpt0==2)stop
!ENDDEBUG

!   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(mkmem==0 .or. squeeze==1)then
     read_cg=0 ; read_cg_disk=1
    end if
#if defined MPI
           if(localrdwf==0)then
            read_cg=0
            read_cg_disk=0

!XG20040106 The following condition is correct
            if((me==0 .and. mpi_enreg%parareel == 0) .or. &
&        (mpi_enreg%me_group_para==0 .and. mpi_enreg%parareel &
&        == 1))read_cg_disk=1
           end if
#endif

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

!DEBUG
!   write(6,*)' wfsinp: before initwf',wff1%offwff
!   write(6,*)' wfsinp: me,read_cg,read_cg_disk=',me,read_cg,read_cg_disk
!   write(6,*)' wfsinp: nban_dp_rdk=',nban_dp_rdk
!ENDDEBUG
    if(read_cg_disk==1)then
     call initwf (cg_disk,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,&
&     isppol0,mcg_disk, &
&     nban_dp_rdk,nkpt0,npw0_k,nspinor0,occ0_k,wff1)
    end if

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

    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(6,*)' wfsinp: me,cg_disk(:,1)=',me,cg_disk(:,1)
!ENDDEBUG

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



#if defined MPI
           if(localrdwf==0)then

          ! Transmit to each of the other processors, when needed
        if (mpi_enreg%parareel == 0) then
            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)

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

               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)

!DEBUG
!               write(6,*)' wfsinp: me,tag,nband_k,npw_k=',me,tag,nband_k,npw_k
!ENDDEBUG

          !     SEND
                if(me==0)then
!BEGIN TF_CHANGES
                   write(6,*)'SENDWFSINP ',me
                 call MPI_SEND(cg_disk,2*npw_k*nspinor*nband_k,&
          &       MPI_DOUBLE_PRECISION,iproc-1,tag,spaceComm,ierr)
!END TF_CHANGES
                end if

          !     RECEIVE
                if(me==iproc-1)then
!BEGIN TF_CHANGES
                 call MPI_RECV(cg_disk,2*npw_k*nspinor*nband_k,&
          &       MPI_DOUBLE_PRECISION,0,tag,spaceComm,status,ierr)
!END TF_CHANGES
                 icg=icg_k(ikpt,isppol)
                 if(squeeze==0)then
                  cg(:,icg+1:icg+npw_k*nspinor*nband_k)=&
          &        cg_disk(:,1:npw_k*nspinor*nband_k)
                 end if
                 ikassoc=ikassoc_me(iproc)
                end if
               end if

              end if
             end do ! iproc
            end if
        else
            if(nproc_max>=2)then
             do iproc=2,nproc_max
          !   Only me=0 and me=iproc-1 are concerned by this
              if(mpi_enreg%me_group_para==0 .or. &
&                 mpi_enreg%me_group_para==iproc-1)then

               ikpt=ikpt_me(iproc)
               isppol=isppol_me(iproc)
!DEBUG
!              write(6,*)' wfsinp: me,ikpt,isppol',me,ikpt,isppol
!ENDDEBUG
               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)
!DEBUG
!               write(6,*)' wfsinp: me,tag,nband_k,npw_k=',me,tag,nband_k,npw_k
!ENDDEBUG
          !     SEND
                if(mpi_enreg%me_group_para==0)then
                 call MPI_SEND(cg_disk,2*npw_k*nspinor*nband_k,&
          & MPI_DOUBLE_PRECISION,iproc-1,tag, &
          & mpi_enreg%kpt_comm_para(mpi_enreg%ipara),ierr)
                end if

          !     RECEIVE
                if(mpi_enreg%me_group_para==iproc-1)then
                 call MPI_RECV(cg_disk,2*npw_k*nspinor*nband_k,&
          & MPI_DOUBLE_PRECISION,0,tag, &
          & mpi_enreg%kpt_comm_para(mpi_enreg%ipara),status,ierr)
                 icg=icg_k(ikpt,isppol)
                 if(squeeze==0)then
                  cg(:,icg+1:icg+npw_k*nspinor*nband_k)=&
          &        cg_disk(:,1:npw_k*nspinor*nband_k)
                 end if
                 ikassoc=ikassoc_me(iproc)
                end if
               end if

              end if
             end do ! iproc
            end if
        end if !parareel

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

          ! Take care of me=0 needing the data
            if((me==0 .and. mpi_enreg%parareel == 0) .or. &
&         (mpi_enreg%me_group_para==0 .and. mpi_enreg%parareel == 1)) 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*nspinor*nband_k+icg)= &
          &     cg_disk(:,1:npw_k*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 !
        if (mpi_enreg%parareel == 0) then
!BEGIN TF_CHANGES
            call MPI_BCAST(eig0_k, &
          &  nban_dp_rdk*(2*nban_dp_rdk)**formeig ,&
          &  MPI_DOUBLE_PRECISION,0,spaceComm,ierr)
!END TF_CHANGES
        else
            call MPI_BCAST(eig0_k, &
          &  nban_dp_rdk*(2*nban_dp_rdk)**formeig,MPI_DOUBLE_PRECISION,&
          &  0,mpi_enreg%kpt_comm_para(mpi_enreg%ipara),ierr)
        end if
           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(6,*)' 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 MPI
           if(localrdwf==0 .and. ikpt/=0)then
           ! Do not forget : ikpt,isppol were redefined ...
            band_index=band_index_k(ikpt,isppol)
!DEBUG
!ENDDEBUG
            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
     write(6,*)' wfsinp: transfer of first-order eigs not yet coded !'
    end if

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

!   Write to disk if appropriate
!   WARNING : mkmem==0 is not allowed yet
!   The coding has to be done ... here, only fragments ...
    if (mkmem==0) then

!    Exclude rwwf time from wfsinp time
     call timab(19,2,tsec)
     if(mpi_enreg%paralbd==0)tim_rwwf=4
     if(mpi_enreg%paralbd==1)tim_rwwf=11
     if(mkmem==0)stop 'not allowed yet'
!    In particular, kg_k is not yet available

     if(master == me ) then
      call rwwf(cg_disk,eig_k,formeig,0,0,ikpt0,isppol0,kg_k,mban_dp_rd,mcg_disk,nband_rdk,&
&      nband_rdk,npw0_k,nspinor0,occ_k,2,1,tim_rwwf,wffnow)
     end if
     call timab(19,1,tsec)
!   Fill the other blocks
    else if(mkmem/=0)then

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

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

!DEBUG
!       write(6,*)' wfsinp: me, for ikassoc,isppol',&
!&       me,ikassoc,isppol
!       write(6,*)' 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 MPI
           if(ikpt_trial/=0)then
             nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
        if (mpi_enreg%parareel == 0) then
             if(minval(abs(mpi_enreg%proc_distrb(ikpt_trial,1:nband_trial,&
&                isppol_trial)-me))/=0)  ikpt_trial=0
        else
             if(mpi_enreg%proc_distrb_para(mpi_enreg%ipara,ikpt_trial) &
&                /= me)  ikpt_trial=0
        end if
           end if
#endif

        if(ikpt_trial/=0 .and. ikpt_trial<=nkpt_eff)then
         write(message,'(a,a,2i5)' )ch10,&
&         ' wfsinp: transfer to ikpt_trial,isppol_trial=',&
&          ikpt_trial,isppol_trial
         call wrtout(6,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*nspinor0*nband_trial+icg_trial)=&
!&          cg(:,1+icg:npw0_k*nspinor0*nband_trial+icg)
          do ii=1,npw0_k*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

!DEBUG
!          if(ikpt_trial==2)then
!           write(6,*)' wfsinp: iband,ipw,cg_disk for ikpt_trial=2'
!           write(6,*)' nban_dp_rdk,npw0_k=', nban_dp_rdk,npw0_k
!           do iband=1,nban_dp_rdk
!            do ipw=1,npw0_k
!             write(6, '(2i5,2es16.6)' )&
!&             iband,ipw,cg_disk(:,ipw+(iband-1)*npw0_k)
!            end do
!           end do
!          end if
!ENDDEBUG


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

!DEBUG
!          write(6,*)' wfsinp: ikpt_trial=',ikpt_trial
!          write(6,*)' nband_k,band_index_trial',nband_k,band_index_trial
!          write(6,*)eig0_k(1:nban_dp_k)
!          write(6,*)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(6,*)' wfsinp: iband,ipw,cg for ikpt_trial=2'
!           write(6,*)' nband_trial,npw_ktrial=',nband_trial,npw_ktrial
!           do iband=1,nband_trial
!            do ipw=1,npw_ktrial
!             write(6, '(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

!DEBUG
!write(6,*)' wfsinp : transferred cg'
!if(ikpt0==2)stop
!ENDDEBUG

!   End if mkmem/=0
    end if
    deallocate(eig_k,eig0_k)
    if(formeig==0)deallocate(occ_k,occ0_k)

!  End the condition of need of this k point
   end if
!DEBUG
!if(me==0 .and. ikpt0==9)then
! write(6,*)' wfsinp : cg array, ikpt0=',ikpt0
! do ipw=1,15
!  write(6, '(i4,2es20.10)' )ipw,cg(:,ipw)
! end do
!end if
!ENDDEBUG


! End of the k loop
  end do
!End of spin loop
 end do
 if(squeeze==1)then
  deallocate(kg_k,kg0_k)
 end if

!DEBUG
!if(me==0)then
! write(6,*)' wfsinp : first band, cg array='
! do ipw=1,15
!  write(6, '(i4,2es20.10)' )ipw,cg(:,ipw)
! end do
!end if
!ENDDEBUG


#if defined MPI
           call timab(67,1,tsec)
        if (mpi_enreg%parareel == 0) then
!BEGIN TF_CHANGES
           call leave_test(mpi_enreg)
!END TF_CHANGES
        end if
           write(message, '(a)' ) ' wfsinp: loop on k-points and spins done in parallel'
           call wrtout(06,message,'COLL')

          !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,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)
        if (mpi_enreg%parareel == 0) then
!BEGIN TF_CHANGES
           call leave_test(mpi_enreg)
!END TF_CHANGES
        end if
           write(message, '(a)' ) ' wfsinp: loop on k-points done in parallel'
           call wrtout(06,message,'COLL')

           call timab(67,2,tsec)

           if(localrdwf==0)then
            deallocate(ikpt_me,nband_k_me)
            deallocate(ikassoc_me,isppol_me)
           end if
#endif

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

!DEBUG
! if(me==0)then
!  write(6,*)' wfsinp : cg array='
!  icg=0
!  do isppol=1,nsppol
!!  do ikpt=1,nkpt
!   do ikpt=1,1
!    nband_k=nband(ikpt+(isppol-1)*nkpt)
!    npw_k=npwarr(ikpt)
!    do iband=1,nband_k
!     write(6,*)' new band, icg=',icg
!     do ipw=1,npw_k
!      write(6, '(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(6,*)' wfsinp : eigen array='
!do ikpt=1,nkpt
! do iband=1,mband
!  write(6,*)'ikpt,iband,eigen',ikpt,iband,eigen(iband+(ikpt-1)*mband)
! end do
!end do
!ENDDEBUG

 if(mkmem/=0)deallocate(icg_k)
 deallocate(band_index_k)

!DEBUG
!write(6,*)' wfsinp : exit'
!write(6,*)cg(:,1)
!write(6,*)cg(:,2)
!stop
!ENDDEBUG

end subroutine wfsinp

!!***

Generated by  Doxygen 1.6.0   Back to index