TABLE OF CONTENTS


ABINIT/inwffil [ Functions ]

[ Top ] [ Functions ]

NAME

 inwffil

FUNCTION

 Do initialization of wavefunction files.
 Also call other relevant routines for this initialisation
 (initialization of wavefunctions from scratch or from file, translations of wavefunctions, ...)

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, AR, MB, MVer)
 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

  ask_accurate= if 1, the wavefunctions and eigenvalues must be
    accurate, that is, they must come from a k point that is
    symmetric of the needed k point, with a very small tolerance,
    the disk file contained sufficient bands to initialize all of them,
    the spinor and spin-polarisation characteristics must be identical
  dtset <type(dataset_type)>=all input variables for this dataset
  ecut=effective kinetic energy planewave cutoff (hartree), beyond
    which the coefficients of plane waves are zero
  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
  hdr <type(hdr_type)>=the header of wf, den and pot files
  ireadwf=option parameter described above for wf initialization
  istwfk(nkpt)=input option parameter that describes the storage of wfs to be initialized here.
  kg(3,mpw*my_nkpt)=dimensionless coords of G vecs in basis sphere at k point
  kptns(3,nkpt)=reduced coords of k points
  localrdwf=(for parallel case) if 1, the wffnm  file is local to each machine
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mkmem=number of k-points in core memory
  mpi_enreg=informations about MPI parallelization
  mpw=maximum number of planewaves as dimensioned in calling routine
  nband(nkpt*nsppol)=number of bands at each k point
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt=number of k points
  npwarr(nkpt)=array holding npw for each k point.
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym=number of symmetry elements in space group
  occ(mband*nkpt*nsppol)=occupations (from disk or left at their initial value)
  optorth= 1 if the WFS have to be orthogonalized; 0 otherwise
  prtvol=control print volume and debugging
  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
  unkg=unit number for storage of basis sphere data: stores indirect
   indexing array and integer coordinates for all planewaves in basis
   sphere for each k point being considered
  unwff1,unwfnow= unit numbers for files wffnm and wft1nm.
  wffnm=name (character data) of file for input wavefunctions.

OUTPUT

  wff1  = structure information for files wffnm .
  wffnow= structure information for wf file wft1nm
  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 (returned if mkmem/=0):
  cg(2,mcg)=complex wf array
    be careful : an array of size cg(2,npw*nspinor), as used
    in the response function code, is not enough !
  wvl <type(wvl_data)>=all wavelets data.

NOTES

 Detailed description :
  Initialize unit wff1%unwff for input of wf data if ireadwf=1
  Opens file on unit wffnow%unwff
   if the storage on disk is needed (mkmem==0)
  Initializes wf data on wffnow%unwff, 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)

 ireadwf options:
   0 => initialize with random numbers or 0 s
   1 => read from disk file wff1, initializing higher bands
        with random numbers or 0 s if not provided in disk file

 The wavefunctions after this initialisation are stored in unit wffnow%unwff

WARNINGS

 The symmetry operations are used to translate the data from one
 k point to another, symmetric, k point.
 They can be completely different from the symmetry operations
 contained on the disk file. No check is performed between the two sets.

 Occupations will not be modified nor output, in the present status of this routine.

 If ground state format (formeig=0) occ(mband*nkpt*nsppol) was output.
 NOT OUTPUT NOW !

PARENTS

      dfpt_looppert,dfptnl_loop,gstate,nonlinear,respfn

CHILDREN

      copy_mpi_enreg,destroy_mpi_enreg,hdr_check,hdr_free,hdr_io,hdr_ncread
      kpgio,listkk,matr3inv,newkpt,pawrhoij_copy,timab,wffopen,wfsinp,wrtout
      wvl_wfsinp_disk,wvl_wfsinp_scratch

SOURCE

116 #if defined HAVE_CONFIG_H
117 #include "config.h"
118 #endif
119 
120 #include "abi_common.h"
121 
122 subroutine inwffil(ask_accurate,cg,dtset,ecut,ecut_eff,eigen,exchn2n3d,&
123 &           formeig,hdr,ireadwf,istwfk,kg,kptns,localrdwf,mband,&
124 &           mcg,mkmem,mpi_enreg,mpw,nband,ngfft,nkpt,npwarr,&
125 &           nsppol,nsym,occ,optorth,symafm,symrel,tnons,unkg,wff1,&
126 &           wffnow,unwff1,wffnm,wvl)
127 
128  use defs_basis
129  use defs_abitypes
130  use defs_wvltypes
131  use m_profiling_abi
132  use m_wffile
133  use m_errors
134  use m_xmpi
135  use m_nctk
136  use m_hdr
137 
138  use m_io_tools, only : file_exists
139  use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_io
140  use m_mpinfo,   only : destroy_mpi_enreg, copy_mpi_enreg
141  use m_kg,       only : kpgio
142 
143 !This section has been created automatically by the script Abilint (TD).
144 !Do not modify the following lines by hand.
145 #undef ABI_FUNC
146 #define ABI_FUNC 'inwffil'
147  use interfaces_14_hidewrite
148  use interfaces_18_timing
149  use interfaces_32_util
150  use interfaces_56_recipspace
151  use interfaces_67_common
152  use interfaces_79_seqpar_mpi, except_this_one => inwffil
153 !End of the abilint section
154 
155  implicit none
156 
157 !Arguments ------------------------------------
158  integer,intent(in) :: ask_accurate,exchn2n3d,formeig,ireadwf,localrdwf,mband,mcg,mkmem,mpw
159  integer,intent(in) :: nkpt,nsppol,nsym,optorth,unkg,unwff1
160  real(dp),intent(in) :: ecut,ecut_eff
161  character(len=*),intent(in) :: wffnm
162  type(MPI_type),intent(inout),target :: mpi_enreg
163  type(dataset_type),intent(in) :: dtset
164  type(hdr_type),intent(inout) :: hdr
165  type(wffile_type),intent(inout) :: wff1
166  type(wffile_type),intent(inout) :: wffnow
167  type(wvl_data),intent(inout) :: wvl
168  integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),ngfft(18)
169  integer,intent(in) :: npwarr(nkpt),symafm(nsym),symrel(3,3,nsym)
170  integer,intent(in),target :: nband(nkpt*nsppol)
171  real(dp),intent(inout),target :: cg(2,mcg),eigen((2*mband)**formeig*mband*nkpt*nsppol)
172  real(dp),intent(in) :: kptns(3,nkpt),tnons(3,nsym)
173  real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
174 
175 !Local variables-------------------------------
176  integer,parameter :: master=0
177  integer :: iomode,accurate,ceksp,debug,doorth,fform,fform_dum,fill
178  integer :: headform0,iband,ibg,ibg0,icg,icg0,icgsft,ieigsft,ierr,ii
179  integer :: ikassoc,ikpt,ikpt0,ikptsp,ikptsp0,imax,increase_nkassoc,isppol,isppol0
180  integer :: mband0,mband0_rd,mband_eff,mcg_disk,me,me0,mkmem0,mpw0
181  integer :: my_nkpt,my_nspinor,my_nspinor0,nband_k,nband0_k
182  integer :: nkassoc,nkpt0,npw,npw0,nspinor0,nspinor_eff,nsppol0,nsppol_eff,nsppol2nspinor
183  integer :: rdwr,randalg,restart,restartpaw,spaceComm,spaceComm_io,sppoldbl,sppoldbl_eff,squeeze
184  logical :: out_of_core
185  real(dp) :: dksqmax,ecut0
186  character(len=500) :: message
187  type(hdr_type) :: hdr0
188  integer :: ngfft0(18)
189  integer,allocatable :: indkk0(:,:),indx(:),istwfk0(:),kg0(:,:)
190  integer,allocatable :: nband0_rd(:),npwarr0(:),npwi(:),npwtot0(:)
191  integer,allocatable,target :: indkk(:,:),nband0(:)
192  integer, pointer :: indkk_eff(:,:),nband_eff(:)
193  logical,allocatable :: my_kpt(:)
194  real(dp) :: gmet(3,3),gmet0(3,3),gprim0(3,3),rprim0(3,3),tsec(2)
195  real(dp),allocatable :: cg_disk(:,:),kptns0(:,:)
196  real(dp),pointer :: cg_eff(:,:),eigen_eff(:)
197  type(MPI_type),pointer :: mpi_enreg0
198 
199 ! *************************************************************************
200 
201  DBG_ENTER("COLL")
202 
203 !Keep track of total time spent in inwffil
204  call timab(710,1,tsec)
205  call timab(711,1,tsec)
206 
207 !Check the validity of formeig
208  if (formeig/=0.and.formeig/=1) then
209    write(message,'(a,i0,a)')' formeig = ',formeig,', but the only allowed values are 0 or 1.'
210    MSG_BUG(message)
211  end if
212 
213 !Init mpi_comm
214  spaceComm=mpi_enreg%comm_cell
215  spaceComm_io=xmpi_comm_self
216  if (mpi_enreg%paral_kgb==1) spaceComm_io= mpi_enreg%comm_bandspinorfft
217  if (mpi_enreg%paral_hf ==1) spaceComm_io= mpi_enreg%comm_hf
218  me=xmpi_comm_rank(spaceComm)
219 
220 !Determine number of k points processed by current node
221  my_nkpt=nkpt;if (size(mpi_enreg%my_kpttab)>0) my_nkpt=maxval(mpi_enreg%my_kpttab)
222  out_of_core=(mkmem==0.and.my_nkpt/=0)
223 
224  ngfft0(:)=ngfft(:)
225  headform0=0 !Default value for headform0 (will be needed later, to read wf blocks)
226 
227 !Chebyshev is more sensitive to the quality of input random numbers, so use a new algorithm
228  if(dtset%wfoptalg == 1) then
229    randalg = 1
230  else
231    ! Otherwise, use compatibility mode
232    randalg = 0
233  end if
234 
235 !If the input data are on disk, determine the kind of restart
236  wff1%fname = wffnm
237 
238 !Checking the existence of data file
239  if (ireadwf==1 .and. .not.file_exists(wff1%fname)) then
240    ! Trick needed to run Abinit test suite in netcdf mode.
241    if (file_exists(nctk_ncify(wff1%fname))) then
242      write(std_out,"(3a)")"- File: ",trim(wff1%fname)," does not exist but found netcdf file with similar name."
243      wff1%fname = nctk_ncify(wff1%fname)
244    end if
245    if (localrdwf/=0 .and. .not. file_exists(wff1%fname)) then
246      MSG_ERROR('Missing data file: '//TRIM(wff1%fname))
247    end if
248  end if
249 
250 !Compute reciprocal space metric gmet
251  call matr3inv(hdr%rprimd,gprim0) ! gprim0 is used as temporary storage
252  gmet=matmul(transpose(gprim0),gprim0)
253 
254  if (ireadwf==1)then
255 
256    iomode=dtset%iomode
257    if (localrdwf==0) then
258      ! This is in case the wff file must be read by only the master proc
259      if (iomode /= IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER
260      !iomode=IO_MODE_FORTRAN_MASTER
261    end if
262 
263    call WffOpen(iomode,spaceComm,wff1%fname,ierr,wff1,master,me,unwff1,spaceComm_io)
264 
265 !  Initialize hdr0 (sent to all procs), thanks to reading of wff1
266    rdwr=1
267    if ( ANY(wff1%iomode == (/IO_MODE_FORTRAN_MASTER, IO_MODE_FORTRAN, IO_MODE_MPI/) )) then
268      call hdr_io(fform_dum,hdr0,rdwr,wff1)
269 #ifdef HAVE_NETCDF
270    else if (wff1%iomode == IO_MODE_ETSF) then
271      call hdr_ncread(hdr0, wff1%unwff, fform_dum)
272 #endif
273    end if
274 
275    ! Handle IO Error.
276    if (fform_dum == 0) then
277      write(message,"(4a)")&
278 &     "hdr_io returned fform == 0 while trying to read the wavefunctions from file: ",trim(wff1%fname),ch10,&
279 &     "This usually means that the file does not exist or that you don't have enough privileges to read it"
280      MSG_ERROR(message)
281    end if
282 
283    call wrtout(std_out,'inwffil: examining the header of disk file '//TRIM(wff1%fname),'COLL')
284 
285 !  Check hdr0 versus hdr (and from now on ignore header consistency and write new info to header for each file)
286    if (dtset%usewvl == 0) then
287 !    wait for plane waves.
288      fform=2
289    else
290 !    wait for wavelets.
291      fform = 200
292    end if
293    call hdr_check(fform,fform_dum,hdr,hdr0,'PERS',restart,restartpaw)
294 
295    nkpt0=hdr0%nkpt
296    nsppol0=hdr0%nsppol
297    headform0=hdr0%headform
298 
299    write(message,'(2a)')'-inwffil : will read wavefunctions from disk file ',trim(wff1%fname)
300    call wrtout(std_out,message,'COLL')
301    call wrtout(ab_out,message,'COLL')
302 
303  else
304    restart=1; restartpaw=0
305 
306 !  Fill some data concerning an hypothetical file to be read
307 !  This is to allow the safe use of same routines than with ireadwf==1.
308    nkpt0=nkpt ; nsppol0=nsppol
309  end if ! end ireadwf
310 
311  sppoldbl=1
312  if(minval(symafm(:))==-1)then
313    if(nsppol0==1 .and. nsppol==2)sppoldbl=2
314  end if
315 
316  ABI_ALLOCATE(indkk,(nkpt*sppoldbl,6))
317  ABI_ALLOCATE(istwfk0,(nkpt0))
318  ABI_ALLOCATE(kptns0,(3,nkpt0))
319  ABI_ALLOCATE(nband0,(nkpt0*nsppol0))
320  ABI_ALLOCATE(npwarr0,(nkpt0))
321 
322  if(restart==2)then ! restart with translations
323 
324    ecut0=hdr0%ecut_eff
325    istwfk0(1:nkpt0)=hdr0%istwfk(1:nkpt0)
326    kptns0(1:3,1:nkpt0)=hdr0%kptns(1:3,1:nkpt0)
327    nband0(1:nkpt0*nsppol0)=hdr0%nband(1:nkpt0*nsppol0)
328    ngfft0(1:3)=hdr0%ngfft(1:3)
329    npwarr0(1:nkpt0)=hdr0%npwarr(1:nkpt0)
330    nspinor0=hdr0%nspinor
331    rprim0(:,:)=hdr0%rprimd(:,:)
332    mpw0=maxval(npwarr0(:))
333 
334 !  Compute reciprocal space metric gmet for unit cell of disk wf
335    call matr3inv(rprim0,gprim0)
336    gmet0=matmul(transpose(gprim0),gprim0)
337 
338    if ((mpi_enreg%paral_kgb==1).or.(mpi_enreg%paral_hf==1)) then
339      ABI_DATATYPE_ALLOCATE(mpi_enreg0,)
340      call copy_mpi_enreg(mpi_enreg,mpi_enreg0)
341      ABI_ALLOCATE(kg0,(3,mpw0*nkpt0))
342      ABI_ALLOCATE(npwtot0,(nkpt0))
343      message="tmpfil"
344      call kpgio(ecut0,dtset%exchn2n3d,gmet0,istwfk0,kg0, &
345 &     kptns0,nkpt0,nband0,nkpt0,'PERS',mpi_enreg0,&
346 &     mpw0,npwarr0,npwtot0,nsppol0)
347 
348      ABI_DEALLOCATE(kg0)
349      ABI_DEALLOCATE(npwtot0)
350    else
351      mpi_enreg0 => mpi_enreg
352    end if
353 
354 !  At this stage, the header of the file wff1i%unwff is read, and
355 !  the pointer is ready to read the first wavefunction block.
356 
357 !  Compute k points from input file closest to the output file
358    call listkk(dksqmax,gmet0,indkk,kptns0,kptns,nkpt0,nkpt,nsym,sppoldbl,symafm,symrel,1)
359 
360  else if (restart==1) then ! direct restart
361 
362 !  Fill variables that must be the same, as determined by hdr_check.f
363 !  This is to allow the safe use of the same routines than with restart==2.
364    nspinor0=dtset%nspinor
365    ecut0=ecut_eff
366    gmet0(:,:)=gmet(:,:)
367    istwfk0(:)=istwfk(:)
368    kptns0(:,:)=kptns(:,:)
369    npwarr0(:)=npwarr(:)
370    mpw0=mpw
371 
372    do isppol=1,sppoldbl
373      do ikpt=1,nkpt
374        indkk(ikpt+(isppol-1)*nkpt,1)=ikpt
375        indkk(ikpt+(isppol-1)*nkpt,2:6)=0
376      end do
377    end do
378    dksqmax=0.0_dp
379 
380 !  The treatment of nband0 asks for some care
381    if(ireadwf==0)then
382      nband0(:)=0
383    else
384      nband0(1:nkpt0*nsppol0)=hdr0%nband(1:nkpt0*nsppol0)
385    end if
386 
387    mpi_enreg0 => mpi_enreg
388 
389  else
390    mpi_enreg0 => mpi_enreg
391  end if
392 
393  if(mpi_enreg0%paral_pert == 1.and.mpi_enreg0%me_pert/=-1) then
394    me0 = mpi_enreg0%me_pert
395  else
396    me0 = mpi_enreg0%me_cell
397  end if
398 
399 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
400 !Before hdr_free:
401 !If restartpaw==1, store hdr0%pawrhoij in hdr%pawrhoij; else if restartpaw==0,
402 !hdr%pawrhoij(:)has been initialized in hdr_init.
403  if(restartpaw==1) then
404    call pawrhoij_copy(hdr0%pawrhoij,hdr%pawrhoij,keep_itypat=.true.)
405  end if
406 
407  call timab(711,2,tsec)
408  call timab(712,1,tsec)
409 
410 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
411 !At this stage, all the relevant information from the header of the disk file,
412 !has been exploited, and stored in variables, on all processors.
413 !It is also contained in hdr0
414 !(on all processors, except if restart=1 and localrdwf=0,
415 !in which case it is only on the master)
416 !These information might be changed later, while processing the
417 !wavefunction data, and converting it. The variable hdr0 might be kept
418 !for further checking, or reference, or debugging, but at present,
419 !it is simpler to close it. The other header, hdr, will be used for the new file, if any.
420 
421  if(ask_accurate==1)then
422 
423 !  Check whether the accuracy requirements might be fulfilled
424    if(ireadwf==0)then
425      write(message,'(9a)')&
426 &     'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
427 &     'present calculation. It was asked that the wavefunctions be accurate,',ch10,&
428 &     'but they were not even read.',ch10,&
429 &     'Action: use a wf file, with ireadwf/=0.'
430      MSG_ERROR(message)
431    end if
432    if(dksqmax>tol12)then
433      write(message, '(9a,es16.6,4a)' )&
434 &     'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
435 &     'present calculation. It was asked that the wavefunctions be accurate, but',ch10,&
436 &     'at least one of the k points could not be generated from a symmetrical one.',ch10,&
437 &     'dksqmax=',dksqmax,ch10,&
438 &     'Action: check your wf file and k point input variables',ch10,&
439 &     '        (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.'
440      MSG_ERROR(message)
441    end if
442    if(dtset%nspinor/=nspinor0)then
443      write(message,'(a,a, a,a,a,a,a, a,a,2i5,a,a)')&
444 &     'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
445 &     'present calculation. It was asked that the wavefunctions be accurate, but',ch10,&
446 &     'nspinor differs in the file from the actual nspinor.',ch10,&
447 &     'nspinor,nspinor0=',dtset%nspinor,nspinor0,ch10,&
448 &     'Action: check your wf file, and nspinor input variables.'
449      MSG_ERROR(message)
450    end if
451    if((nsppol>nsppol0 .and. sppoldbl==1) .or. nsppol<nsppol0 ) then
452      write(message,'(a,a, a,a,a,a,a, a,a,3i5,a,a)')&
453 &     'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
454 &     'present calculation. It was asked that the wavefunctions be accurate, but',ch10,&
455 &     'the nsppol variables do not match in the file and in the actual calculation',ch10,&
456 &     'nsppol,nsppol,sppoldbl=',dtset%nspinor,nspinor0,sppoldbl,ch10,&
457 &     'Action: check your wf file, and nsppol input variables.'
458      MSG_ERROR(message)
459    end if
460 
461 !  Now, check the number of bands
462    accurate=1
463    do isppol=1,nsppol
464      do ikpt=1,nkpt
465        ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1)
466        ikptsp =ikpt +(isppol-1)*nkpt
467        ikptsp0=ikpt0+(isppol-1)*(2-sppoldbl)*nkpt0
468        if(nband0(ikptsp0)<nband(ikptsp))accurate=0
469      end do
470    end do
471    if(accurate==0)then
472      write(message,'(a,a, a,a,a,a,a, a,a)')&
473 &     'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
474 &     'present calculation. It was asked that the wavefunctions be accurate,',ch10,&
475 &     'but the number of bands differ in the file and in the actual calculation.',ch10,&
476 &     'Action: use a wf file with the correct characteristics.'
477      MSG_ERROR(message)
478    end if
479 
480  end if
481 
482 !Flag: do we need to translate WF to (from) spinors ?
483  nsppol2nspinor=0
484  if (nsppol0==2.and.dtset%nspinor==2) nsppol2nspinor=+1
485  if (nspinor0==2.and.nsppol==2) nsppol2nspinor=-1
486 
487 !Take into account parallism over spinors
488  my_nspinor =max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
489  my_nspinor0=max(1,nspinor0/mpi_enreg0%nproc_spinor)
490 
491 !Not all bands might be read, if not needed to fill the wavefunctions
492  mband0=maxval(nband0(1:nkpt0*nsppol0))
493  mband0_rd=min(mband0,(mband/dtset%nspinor)*nspinor0)
494 
495 !****************************************************************************
496 !If needed, transfer the input wf from disk to core memory
497 !(in the parallel case, it allows to change localrdwf=0 in localrdwf=1)
498 
499  mkmem0=0
500 
501  if(xmpi_paral == 1 .or. mpi_enreg%paral_kgb == 1 .or. mpi_enreg%paral_hf == 1) then
502    if(localrdwf==0 .and. out_of_core)then
503      MSG_BUG('localrdwf==0 and mkmem==0 (out-of-core solution) are not allowed together (yet)')
504    end if
505  end if
506 
507  call timab(712,2,tsec)
508 
509 !Here, treat reading wavefunctions with mkmem/=0, first step
510  if(ireadwf==1 .and. (.not.out_of_core))then
511 
512    call timab(713,1,tsec)
513 
514 !  if(restart==1 .and. ireadwf==1 .and. mkmem/=0)then
515 
516 !  Compute table of k point associations. Make a trial choice for nkassoc.
517    nkassoc=(nkpt/nkpt0+1)*2
518    ABI_ALLOCATE(indkk0,(nkpt0,nkassoc))
519 !  Infinite loops are allowed in F90
520    do
521      indkk0(:,:)=0
522      increase_nkassoc=0
523      do ikpt=1,nkpt*sppoldbl
524        ikpt0=indkk(ikpt,1)
525        do ikassoc=1,nkassoc
526          if(indkk0(ikpt0,ikassoc)==0)then
527            indkk0(ikpt0,ikassoc)=ikpt
528            exit
529          end if
530          if(nkassoc==ikassoc)increase_nkassoc=1
531        end do
532        if(increase_nkassoc==1)then
533          ABI_DEALLOCATE(indkk0)
534          nkassoc=2*nkassoc
535          ABI_ALLOCATE(indkk0,(nkpt0,nkassoc))
536          exit
537        end if
538      end do
539      if(increase_nkassoc==0)exit
540    end do
541 
542 !  DEBUG
543 !  write(std_out,*)' inwffil: indkk0, nkassoc=',nkassoc
544 !  do ikpt0=1,nkpt0
545 !  write(std_out,*)' ikpt0,indkk0(ikpt0,1)=',ikpt0,indkk0(ikpt0,1)
546 !  end do
547 !  ENDDEBUG
548 
549 !  DEBUG
550 !  write(std_out,*)' inwffil : indkk(:,1)=',indkk(:,1)
551 !  write(std_out,*)' inwffil : sppoldbl=',sppoldbl
552 !  ENDDEBUG
553 
554 !  To treat the case (nsppol0=2,nspinor0=1)<->(nsppol=1,nspinor=2),
555 !  apply the following trick:
556 !  1- We call wfsinp with fake arguments (nsppol_eff and nspinor_eff)
557 !  2- We transform collinear polarized WF into spinors
558 !  or  spinors into collinear polarized WF
559    if (nsppol2nspinor/=0.and.out_of_core.and.dtset%usewvl==0) then
560      write(message, '(7a)')&
561 &     'When mkmem=0 (out-of-core), the wavefunction translator is unable',ch10,&
562 &     'to interchange spin-polarized wfs and spinor wfs.',ch10,&
563 &     'Action: use a non-spin-polarized wf to start a spinor wf,',ch10,&
564 &     '        and a non-spinor wf to start a spin-polarized wf.'
565      MSG_ERROR(message)
566    end if
567 
568 !  === Fake arguments definition for wfsinp
569    if (nsppol2nspinor==0.or.dtset%usewvl/=0) then
570      indkk_eff => indkk
571      nband_eff => nband
572      eigen_eff => eigen
573      cg_eff => cg
574      nspinor_eff=dtset%nspinor;nsppol_eff=nsppol;sppoldbl_eff=sppoldbl
575      mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff))
576    else if (nsppol2nspinor==1.and.(.not.out_of_core)) then
577      nsppol_eff=2;nspinor_eff=1;sppoldbl_eff=1
578      ABI_ALLOCATE(indkk_eff,(nkpt*sppoldbl_eff,6))
579      ABI_ALLOCATE(nband_eff,(nkpt*nsppol_eff))
580      indkk_eff(1:nkpt,1:6)   =indkk(1:nkpt,1:6)
581      nband_eff(1:nkpt)       =nband(1:nkpt)/2
582      nband_eff(1+nkpt:2*nkpt)=nband(1:nkpt)/2
583      mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff))
584      eigen_eff => eigen
585      cg_eff => cg
586    else if (nsppol2nspinor==-1.and.(.not.out_of_core)) then
587 !    WARNING: MT 07072011 -> this is memory consuming
588 !    A copy a spinorial WF (and eigenvalues) is temporary kept in memory;
589 !    But the case (nspinor=2 => nsppol=2) might be rare
590 !    and only useful for testing purposes.
591 !    => print a warning for the user
592 !    NOTE: in that case (nsppol=2), parallelization over spinors is not activated
593 
594      write(message,'(5a)')&
595 &     'In the case of spinor WF read from disk and converted into',ch10,&
596 &     'spin-polarized non-spinor WF, the WF translator is memory',ch10,&
597 &     'consuming (a copy a spinor WF is temporary stored in memory).'
598      MSG_WARNING(message)
599 
600      nsppol_eff=1;nspinor_eff=2;sppoldbl_eff=1
601      ABI_ALLOCATE(indkk_eff,(nkpt*sppoldbl_eff,6))
602      ABI_ALLOCATE(nband_eff,(nkpt*nsppol_eff))
603      indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6)
604      nband_eff(1:nkpt)    =2*nband(1:nkpt)
605      mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff))
606      ABI_ALLOCATE(eigen_eff,((2*mband_eff)**formeig*mband_eff*nkpt*nsppol_eff))
607      ABI_ALLOCATE(cg_eff,(2,mpw0*nspinor_eff*mband_eff*mkmem*nsppol_eff))
608    end if
609 
610 !  === nband0 argument definition for wfsinp
611    squeeze=0
612    ABI_ALLOCATE(cg_disk,(0,0))
613    if(.not.out_of_core)then
614      ABI_ALLOCATE(nband0_rd,(nkpt0*nsppol0))
615      nband0_rd(:)=0
616      do isppol=1,nsppol_eff
617        do ikpt=1,nkpt
618          ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1)
619          isppol0=min(isppol,nsppol0)
620          ikptsp =ikpt +(isppol -1)*nkpt
621          ikptsp0=ikpt0+(isppol0-1)*(2-sppoldbl)*nkpt0
622          nband0_k=min(nband0(ikptsp0),(nband_eff(ikptsp)/nspinor_eff)*nspinor0)
623          nband0_rd(ikptsp0)=max(nband0_rd(ikptsp0),nband0_k)
624          npw0=npwarr0(ikpt0)
625          npw =npwarr (ikpt)
626          if(npw0*nspinor0*nband0_k > npw*nspinor_eff*nband_eff(ikptsp))squeeze=1
627        end do
628      end do
629      if(squeeze==1)then
630        mcg_disk=mpw0*my_nspinor0*mband0_rd
631        ABI_DEALLOCATE(cg_disk)
632        ABI_ALLOCATE(cg_disk,(2,mcg_disk))
633      else
634        if(xmpi_paral == 1 .or. mpi_enreg0%paral_kgb == 1 .or. mpi_enreg0%paral_hf == 1)then
635          if(localrdwf==0)then
636            mcg_disk=mpw0*my_nspinor0*mband0_rd
637            ABI_DEALLOCATE(cg_disk)
638            ABI_ALLOCATE(cg_disk,(2,mcg_disk))
639          end if
640        end if
641      end if
642    end if
643 
644    call timab(713,2,tsec)
645    call timab(714,1,tsec)
646 
647 !  === call to wfsinp
648    if (dtset%usewvl == 0) then
649      call wfsinp(cg_eff,cg_disk,ecut,ecut0,ecut_eff,eigen,&
650 &     exchn2n3d,formeig,gmet,gmet0,headform0,&
651 &     indkk_eff,indkk0,istwfk,istwfk0,kptns,kptns0,localrdwf,&
652 &     mband_eff,mcg,mcg_disk,mpi_enreg,mpi_enreg0,mpw,mpw0,&
653 &     nband_eff,nband0_rd,ngfft,nkassoc,nkpt,nkpt0,npwarr,npwarr0,nspinor_eff,nspinor0,&
654 &     nsppol_eff,nsppol0,nsym,occ,optorth,dtset%prtvol,randalg,restart,hdr%rprimd,sppoldbl_eff,squeeze,&
655 &     symrel,tnons,wff1)
656      if (nsppol2nspinor/=0)  then
657        ABI_DEALLOCATE(indkk_eff)
658        ABI_DEALLOCATE(nband_eff)
659      end if
660    else
661 !    Read wavefunctions from file.
662      call wvl_wfsinp_disk(dtset, hdr0, hdr, mpi_enreg, occ, 1, &
663 &     hdr%rprimd, wff1, wvl%wfs, wvl%descr, hdr%xred)
664    end if
665 
666    call timab(714,2,tsec)
667    call timab(715,1,tsec)
668 
669 !  Now, update xyz0 variables, for use in newkpt
670    nband0(:)=nband0_rd(:)
671 
672 !  If squeeze, the conversion was done in wfsinp, so no conversion left.
673    if(squeeze==1)then
674      ecut0=ecut_eff
675      gmet0(:,:)=gmet(:,:)
676      ABI_DEALLOCATE(kptns0)
677      ABI_DEALLOCATE(istwfk0)
678      ABI_DEALLOCATE(nband0)
679      ABI_DEALLOCATE(npwarr0)
680      ABI_ALLOCATE(kptns0,(3,nkpt))
681      ABI_ALLOCATE(istwfk0,(nkpt))
682      ABI_ALLOCATE(nband0,(nkpt*nsppol))
683      ABI_ALLOCATE(npwarr0,(nkpt))
684      kptns0(:,:)=kptns(:,:)
685      istwfk0(:)=istwfk(:)
686      npwarr0(:)=npwarr(:)
687      nband0(:)=0
688      do isppol=1,nsppol
689        do ikpt=1,nkpt
690          ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1)
691          isppol0=min(isppol,nsppol0)
692          ikptsp =ikpt +(isppol -1)*nkpt
693          ikptsp0=ikpt0+(isppol0-1)*(sppoldbl-1)*nkpt0
694          nband0(ikptsp)=(nband0_rd(ikptsp0)/nspinor0)*dtset%nspinor
695        end do
696      end do
697      do ikpt=1,nkpt
698        indkk(ikpt,1)=ikpt
699        indkk(ikpt,2:6)=0
700      end do
701 !    This transfer must come after the nband0 transfer
702      nspinor0=dtset%nspinor
703      nkpt0=nkpt
704      nsppol0=nsppol
705    end if ! end squeeze == 1
706 
707 !  The input wavefunctions have been transferred from disk to core memory
708    mkmem0=mkmem
709 
710    ABI_DEALLOCATE(indkk0)
711    ABI_DEALLOCATE(nband0_rd)
712    ABI_DEALLOCATE(cg_disk)
713 
714    call timab(715,2,tsec)
715 
716  else !ireadwf == 0
717    if (dtset%usewvl == 1) then
718 
719      call timab(714,1,tsec)
720 !    Compute wavefunctions from input guess.
721      call wvl_wfsinp_scratch(dtset, mpi_enreg, occ, hdr%rprimd, wvl, hdr%xred)
722      call timab(714,2,tsec)
723    end if
724  end if
725 
726  call timab(716,1,tsec)
727 
728 !=== Eventual conversion of WF into (from) spinors
729  if (dtset%usewvl==0) then
730 
731 !  ***** No conversion (standard case) ****
732    if (nsppol2nspinor==0) then
733      nspinor_eff=nspinor0;nsppol_eff=nsppol0;sppoldbl_eff=sppoldbl
734      indkk_eff => indkk
735      nband_eff => nband0
736 
737 !    ***** Conversion from collinear to spinorial WF ****
738    else if (nsppol2nspinor==1.and.(.not.out_of_core)) then
739 !    Translate the WF and eigenvalues from nsppol=2 to nspinor=2
740 !    This is tricky (because we do not want to create a temporary array for cg)
741      nsppol_eff=1;nspinor_eff=2;sppoldbl_eff=1
742      ABI_ALLOCATE(indkk_eff,(nkpt*sppoldbl_eff,6))
743      ABI_ALLOCATE(nband_eff,(nkpt0*nsppol_eff))
744      indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6)
745      nband_eff(1:nkpt0)=2*nband0(1:nkpt0)
746 !    Compute some shifts from isspol0=1 to isppol0=2
747      imax=0;icgsft=0;ieigsft=0
748      ABI_ALLOCATE(my_kpt,(nkpt0))
749      do ikpt0=1,nkpt0
750        nband0_k=nband0(ikpt0);nband_k=nband(ikpt0)
751        my_kpt(ikpt0)=(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me0)))
752        ieigsft=ieigsft+(2*nband0_k)**formeig*nband0_k
753        if(my_kpt(ikpt0)) then
754          imax=imax+nband0_k;icgsft=icgsft+nband0_k*npwarr0(ikpt0)
755        end if
756      end do
757 !    --- First version: no parallelization over spinors
758      if (mpi_enreg0%paral_spinor==0) then
759 !      Compute some useful indexes
760        ABI_ALLOCATE(indx,(2*imax))
761        ABI_ALLOCATE(npwi,(imax))
762        ii=0;icg=0
763        do ikpt0=1,nkpt0
764          if(my_kpt(ikpt0)) then
765            nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
766            do iband=1,nband0_k
767              ii=ii+1;npwi(ii)=npw0
768              indx(2*ii-1)=icg+mpw0;indx(2*ii)=icg+2*mpw0
769              icg=icg+4*mpw0
770            end do
771          end if
772        end do
773 !      Expand WF in cg (try to use the whole array)
774        ii=nsppol0*imax;icg0=nsppol0*icgsft
775        do isppol=nsppol0,1,-1
776          do ikpt0=nkpt0,1,-1
777            if(my_kpt(ikpt0)) then
778              nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
779              do iband=nband0_k,1,-1
780                icg0=icg0-npw0
781                if (indx(ii)<icg0) then
782                  MSG_BUG("Unable to read WF!")
783                end if
784                cg(:,indx(ii)+1:indx(ii)+npw0)=cg(:,icg0+1:icg0+npw0)
785                ii=ii-1
786              end do
787            end if
788          end do
789        end do
790 !      Convert polarized WF into spinors
791        ii=1
792        do ikpt0=1,nkpt0
793          if(my_kpt(ikpt0)) then
794            nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
795            do iband=1,nband0_k
796              npw0=npwi(ii)
797              cg(:,indx(2*ii-1)-mpw0+1:indx(2*ii-1)-mpw0+npw0)=cg(:,indx(ii)+1:indx(ii)+npw0)
798              cg(:,indx(2*ii  )+mpw0+1:indx(2*ii  )+mpw0+npw0)=cg(:,indx(ii+imax)+1:indx(ii+imax)+npw0)
799              ii=ii+1
800            end do
801          end if
802        end do
803 !      Compress new cg array (from mpw to npw) and cancel zero-components
804        icg0=0;icg=0
805        do ikpt0=1,nkpt0
806          if(my_kpt(ikpt0)) then
807            nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
808            do iband=1,nband0_k
809              cg(:,icg0       +1:icg0+  npw0)=cg(:,icg+1:icg+npw0)
810              cg(:,icg0+  npw0+1:icg0+2*npw0)=zero
811              cg(:,icg0+2*npw0+1:icg0+3*npw0)=zero
812              cg(:,icg0+3*npw0+1:icg0+4*npw0)=cg(:,icg+3*mpw0+1:icg+3*mpw0+npw0)
813              icg0=icg0+4*npw0;icg=icg+4*mpw0
814            end do
815          end if
816        end do
817 !      --- Second version: parallelization over spinors
818      else
819 !      Compute some useful indexes
820        ABI_ALLOCATE(indx,(imax))
821        ABI_ALLOCATE(npwi,(imax))
822        ii=0;icg=0
823        do ikpt0=1,nkpt0
824          if(my_kpt(ikpt0)) then
825            nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
826            do iband=1,nband0_k
827              ii=ii+1;npwi(ii)=npw0
828              indx(ii)=icg+mpi_enreg0%me_spinor*mpw0
829              icg=icg+2*mpw0
830            end do
831          end if
832        end do
833 !      Expand WF in cg
834        ii=(mpi_enreg0%me_spinor+1)*imax;icg0=(mpi_enreg0%me_spinor+1)*icgsft
835        do ikpt0=nkpt0,1,-1
836          if(my_kpt(ikpt0)) then
837            nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
838            do iband=nband0_k,1,-1
839              icg0=icg0-npw0
840              if (indx(ii)<icg0) then
841                MSG_BUG("Unable to read WF!")
842              end if
843              cg(:,indx(ii)+1:indx(ii)+npw0)=cg(:,icg0+1:icg0+npw0)
844              ii=ii-1
845            end do
846          end if
847        end do
848 !      Compress new cg array (from mpw to npw) and cancel zero-components
849        icg0=0;icg=0
850        do ikpt0=1,nkpt0
851          if(my_kpt(ikpt0)) then
852            nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
853            do iband=1,nband0_k
854              if (mpi_enreg0%me_spinor==0) then
855                cg(:,icg0     +1:icg0+  npw0)=cg(:,icg+1:icg+npw0)
856                cg(:,icg0+npw0+1:icg0+2*npw0)=zero
857              else
858                cg(:,icg0     +1:icg0+  npw0)=zero
859                cg(:,icg0+npw0+1:icg0+2*npw0)=cg(:,icg+mpw0+1:icg+mpw0+npw0)
860              end if
861              icg0=icg0+2*npw0;icg=icg+2*mpw0
862            end do
863          end if
864        end do
865      end if
866 !    Translate eigenvalues
867      ibg0=2*ieigsft;ibg=2*ieigsft
868      do ikpt0=nkpt0,1,-1
869        nband0_k=nband0(ikpt0)
870        ibg0=ibg0-  nband0_k*(2*nband0_k)**formeig
871        ibg =ibg -2*nband0_k*(2*nband0_k)**formeig
872        if(my_kpt(ikpt0)) then
873          do iband=nband0_k*(2*nband0_k)**formeig,1,-1
874            eigen(2*iband-1+ibg)=eigen(iband+ibg0-ieigsft)
875            eigen(2*iband  +ibg)=eigen(iband+ibg0)
876          end do
877        end if
878      end do
879      ABI_DEALLOCATE(indx)
880      ABI_DEALLOCATE(npwi)
881      ABI_DEALLOCATE(my_kpt)
882 
883 !    ***** Conversion from spinorial to collinear WF ****
884    else if (nsppol2nspinor==-1.and.(.not.out_of_core)) then
885 !    In that case parallelization over spinors is never activated
886      nsppol_eff=2;nspinor_eff=1;sppoldbl_eff=1
887      ABI_ALLOCATE(indkk_eff,(nkpt*sppoldbl_eff,6))
888      ABI_ALLOCATE(nband_eff,(nkpt0*nsppol_eff))
889      indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6)
890      nband_eff(1:nkpt0)        =nband0(1:nkpt0)/2
891      nband_eff(1+nkpt0:2*nkpt0)=nband0(1:nkpt0)/2
892 !    Compute shifts from isspol0=1 to isppol0=2
893      icgsft=0;ieigsft=0
894      do ikpt0=1,nkpt0
895        nband0_k=nband0(ikpt0);nband_k=nband(ikpt0)
896        ieigsft=ieigsft+(nband0_k/2)*(nband0_k)**formeig
897        if(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me))) &
898 &       icgsft=icgsft+(nband0_k/2)*npwarr0(ikpt0)
899      end do
900 !    Translate the WF and eigenvalues from nspinor=2 to nsppol=2
901      icg0=0;icg=0;ibg=0
902      do ikpt0=1,nkpt0
903        nband0_k=nband0(ikpt0);nband_k=nband(ikpt0);npw0=npwarr0(ikpt0)
904        if(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me))) then
905          do iband=1,nband0_k/2
906            do ii=1,npw0
907              cg(:,ii+icg)       =cg_eff(:,ii+icg0)
908              cg(:,ii+icg+icgsft)=cg_eff(:,ii+icg0+3*npw0)
909            end do
910            icg0=icg0+4*npw0;icg=icg+npw0
911          end do
912          do iband=(nband0_k/2)*(nband0_k)**formeig,1,-1
913            eigen(iband+ibg)        =eigen_eff(2*iband-1+2*ibg)
914            eigen(iband+ibg+ieigsft)=eigen_eff(2*iband  +2*ibg)
915 !          occ(iband+ibg)        =occ_eff(2*iband-1+2*ibg)
916 !          occ(iband+ibg+ieigsft)=occ_eff(2*iband  +2*ibg)
917          end do
918        end if
919        ibg=ibg+(nband0_k/2)*(nband0_k)**formeig
920      end do
921      ABI_DEALLOCATE(cg_eff)
922      ABI_DEALLOCATE(eigen_eff)
923 
924    else
925      MSG_BUG('unable to interchange nsppol and nspinor when mkmem=0')
926    end if
927  end if
928 
929 !Clean hdr0
930  !if (ireadwf==1)then
931  !  if( restart==2 .or. localrdwf==1 .or. master==me)then
932  !    call hdr_free(hdr0)
933  !  end if
934  !end if
935  call hdr_free(hdr0)
936 
937  call timab(716,2,tsec)
938  call timab(717,1,tsec)
939 
940 
941 !****************************************************************************
942 !Now, treat translation of wavefunctions if wavefunctions are planewaves
943 
944  ceksp=0; debug=0; doorth=1; fill=1
945  if (dtset%usewvl == 0) then
946 
947    call newkpt(ceksp,cg,debug,ecut0,ecut,ecut_eff,eigen,exchn2n3d,&
948 &   fill,formeig,gmet0,gmet,headform0,indkk_eff,&
949 &   ab_out,ireadwf,istwfk0,istwfk,kg,kptns0,kptns,&
950 &   mband,mcg,mkmem0,mkmem,mpi_enreg0,mpi_enreg,&
951 &   mpw0,mpw,my_nkpt,nband_eff,nband,ngfft0,ngfft,nkpt0,nkpt,npwarr0,npwarr,&
952 &   nspinor_eff,dtset%nspinor,nsppol_eff,nsppol,nsym,occ,optorth,&
953 &   dtset%prtvol,randalg,restart,hdr%rprimd,sppoldbl_eff,symrel,tnons,unkg,wff1,wffnow)
954 
955    if (nsppol2nspinor/=0)  then
956      ABI_DEALLOCATE(nband_eff)
957    end if
958 
959  end if ! dtset%usewvl == 0
960 
961 !****************************************************************************
962 
963  ABI_DEALLOCATE(indkk)
964  ABI_DEALLOCATE(istwfk0)
965  ABI_DEALLOCATE(kptns0)
966  ABI_DEALLOCATE(nband0)
967  ABI_DEALLOCATE(npwarr0)
968  if (restart==2 .and.(mpi_enreg0%paral_kgb==1 .or. mpi_enreg0%paral_hf == 1)) then
969    call destroy_mpi_enreg(mpi_enreg0)
970    ABI_DATATYPE_DEALLOCATE(mpi_enreg0)
971  else
972    nullify(mpi_enreg0)
973  end if
974 
975  call timab(717,2,tsec)
976  call timab(710,2,tsec)
977 
978  DBG_EXIT("COLL")
979 
980 end subroutine inwffil