TABLE OF CONTENTS


ABINIT/newkpt [ Functions ]

[ Top ] [ Functions ]

NAME

 newkpt

FUNCTION

 This subroutine writes a starting guess for wave function (set 2)
 It performs a "zero order" interpolation, ie simply
 searches the nearest available k-point.
 The data (set 1) associated with this point is either
 read from a disk file (with a random access reading routine),
 or input as argument.

COPYRIGHT

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

  ceksp2=if 1, center the sphere of pw on Gamma; if 0, on each k-point.
  doorth=1 to do orthogonalization
  debug=>0 for debugging output
  ecut1=kinetic energy cutoffs for basis sphere 1 (hartree)
  ecut2=kinetic energy cutoffs beyond which the coefficients of wf2 vanish (Ha)
  ecut2_eff=kinetic energy cut-off for basis sphere 2 (hartree)
  exchn2n3d=if 1, n2 and n3 are exchanged
  fill=if 1, fill the supplementary bands ; if 0, reduce the number of bands
             Note : must have fill/=0 in the parallel execution
  formeig=if 0, GS format for wfs, eig and occ ; if 1, RF format.
  gmet1(3,3), gmet2(3,3)=reciprocal space metrics (bohr^-2)
  headform1=header format (might be needed to read the block of wfs)
  indkk(nkpt2*sppoldbl,6)=describe k point number of kptns1 that allows to
   generate wavefunctions closest to given kpt2 (and possibly isppol2=2)
   indkk(:,1)=k point number of kpt1
   indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a
    (if 0, means no symmetry operation, equivalent to identity )
   indkk(:,3:5)=shift in reciprocal space to be given to kpt1a,
    to give kpt1b, that is the closest to ikpt2.
   indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
  iout=unit number for output file
  ireadwf=if 0, no reading of disk wavefunction file (random or 0.0 initialisation)
  istwfk1(nkpt1)=input parameter that describes the storage of wfs in set1
  istwfk2(nkpt2)=input parameter that describes the storage of wfs in set2
  kg2(3,mpw2*mkmem2)=dimensionless coords of G vecs in basis sphere at k point
  kptns1(3,nkpt1), kptns2(3,nkpt2)=k point sets (reduced coordinates)
  mband2= maximum number of bands of the output wavefunctions
  mcg=dimension of the cg array
   In case mkmem2/=0, all the output data must find their place in cg,
    so that mcg must be at least Sum(ikpt,isppol) [npw*nspinor*nband](ikpt,isppol)
    where these data are related to the output parameters
   In case mkmem1/=0, the same is true, for the input parameters,
    however, the maximum number of bands that will be read
    will be at most (mband2/nspinor2)*nspinor1
   In case mkmem1==0 and mkmem2==0, one must have at least mpw*nspinor*mband
    for BOTH the input and output parameters, taking into account the
    maximal number of band to be read, described above.
   In case mkmem1/=0 and mkmem2/=0, it is expected that the input cg array
    is organised using the output parameters nkpt2, nband2 ...
    This is needed, in order to use the same pointer.
  mkmem1= if 0, the input wf, eig, occ are available from disk
  mkmem2= if 0, the output wf, eig, occ must be written onto disk
  mpi_enreg1=informations about MPI parallelization, for the input wf file
  mpi_enreg2=informations about MPI parallelization, for the output wf file
  mpw1=maximum allowed number of planewaves at any k, for the input wf file
  mpw2=maximum allowed number of planewaves at any k, for the output wf file
  my_nkpt2= number of k points for the output wf file, handled by current processus
  nband1(nkpt1*nsppol1)=number of bands, at each k point, on disk
  nband2(nkpt2*nsppol2)=desired number of bands at each k point
  ngfft1(18)=all needed information about 3D FFT, for the input wf file
  ngfft2(18)=all needed information about 3D FFT, for the output wf file
             see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt1, nkpt2=number of k points in each set
  npwarr1(nkpt1)=array holding npw for each k point (input wf file).
  npwarr2(nkpt2)=array holding npw for each k point (output wf file).
  nspinor1,nspinor2=number of spinorial components of the wavefunctions
   for each wf file (input or output)
  nsppol1=1 for unpolarized, 2 for spin-polarized, input wf file
  nsppol2=1 for unpolarized, 2 for spin-polarized, output wf file
  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
  randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility
  restart= if 2, 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
  symrel(3,3,nsym)=symmetry operations in real space in terms
   of primitive translations
  tnons(3,nsym)=nonsymmorphic translations for symmetry operations
  unkg2=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 (kptns2 set)
  wffinp=structure info of input wf file unit number
  wffout=structure info of output wf file unit number

OUTPUT

  (see side effects)

SIDE EFFECTS

     The following arrays are input if mkmem1/=0, otherwise their input
     values are taken from disk, and are output if mkmem2/=0, otherwise
     their output values are written on disk.
     The location of the block for a given spin-k point at input MUST
     be the same as the location of the corresponding spin-k point at output.
  cg(2,mcg)=complex wf array
  eigen(mband2*(2*mband2)**formeig *nkpt2*nsppol2)=
    eigenvalues (input or init to large number for GS or init to 0.0 for RF), (Ha)
  occ(mband2*nkpt2*nsppol2)=occupation (input or init to 0.0)  NOT USED NOW

NOTES

 * When reading from disk, it is expected that the next record of
 the wffinp%unwff disk unit is the first record of the first wavefunction block.

 * When the data is input as argument, it is assumed that the
 data for each spin- k wavefunction block is located at the proper
 corresponding location of the output array (this is to be described).

 * The information is pumped onto an fft box for the conversion.
 This allows for changing the number of plane waves.

 * In the present status of this routine, occ is not output.

PARENTS

      inwffil

CHILDREN

      pareigocc,prmat,randac,rdnpw,rwwf,timab,wfconv,wffreadskipk,wrtout

SOURCE

134 #if defined HAVE_CONFIG_H
135 #include "config.h"
136 #endif
137 
138 #include "abi_common.h"
139 
140 subroutine newkpt(ceksp2,cg,debug,ecut1,ecut2,ecut2_eff,eigen,exchn2n3d,fill,&
141 &                  formeig,gmet1,gmet2,headform1,indkk,iout,ireadwf,&
142 &                  istwfk1,istwfk2,kg2,kptns1,kptns2,mband2,mcg,mkmem1,mkmem2,&
143 &                  mpi_enreg1,mpi_enreg2,mpw1,mpw2,my_nkpt2,nband1,nband2,&
144 &                  ngfft1,ngfft2,nkpt1,nkpt2,npwarr1,npwarr2,nspinor1,nspinor2,&
145 &                  nsppol1,nsppol2,nsym,occ,optorth,prtvol,randalg,restart,rprimd,&
146 &                  sppoldbl,symrel,tnons,unkg2,wffinp,wffout)
147 
148  use defs_basis
149  use defs_abitypes
150  use m_errors
151  use m_profiling_abi
152  use m_wffile
153  use m_xmpi
154 
155  use m_pptools,    only : prmat
156  use m_occ,        only : pareigocc
157 
158 !This section has been created automatically by the script Abilint (TD).
159 !Do not modify the following lines by hand.
160 #undef ABI_FUNC
161 #define ABI_FUNC 'newkpt'
162  use interfaces_14_hidewrite
163  use interfaces_18_timing
164  use interfaces_32_util
165  use interfaces_56_io_mpi
166  use interfaces_62_iowfdenpot
167  use interfaces_66_wfs
168 !End of the abilint section
169 
170  implicit none
171 
172 !Arguments ------------------------------------
173 !scalars
174  integer,intent(in) :: ceksp2,debug,exchn2n3d,fill,formeig,headform1,iout
175  integer,intent(in) :: ireadwf,mband2,mcg,mkmem1,mkmem2,mpw1,mpw2,my_nkpt2,nkpt1,nkpt2
176  integer,intent(in) :: nspinor1,nspinor2,nsppol1,nsppol2,nsym,optorth,prtvol,restart
177  integer,intent(in) :: randalg,sppoldbl,unkg2
178  real(dp),intent(in) :: ecut1,ecut2,ecut2_eff
179  type(MPI_type),intent(inout) :: mpi_enreg1,mpi_enreg2
180  type(wffile_type),intent(inout) :: wffinp,wffout
181 !arrays
182  integer,intent(in) :: indkk(nkpt2*sppoldbl,6),istwfk1(nkpt1),istwfk2(nkpt2)
183  integer,intent(in) :: kg2(3,mpw2*mkmem2),nband1(nkpt1*nsppol1)
184  integer,intent(in) :: nband2(nkpt2*nsppol2),ngfft1(18),ngfft2(18)
185  integer,intent(in) :: npwarr1(nkpt1),npwarr2(nkpt2),symrel(3,3,nsym)
186  real(dp),intent(in) :: gmet1(3,3),gmet2(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2)
187  real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym)
188  real(dp),intent(inout) :: cg(2,mcg) !vz_i pw_orthon vecnm
189  real(dp),intent(inout) :: eigen(mband2*(2*mband2)**formeig*nkpt2*nsppol2)!vz_i newocc
190  real(dp),intent(inout) :: occ(mband2*nkpt2*nsppol2) !vz_i
191 
192 !Local variables-------------------------------
193 !scalars
194  integer,parameter :: init_random=-5,nkpt_max=50,tobox=1,tosph=-1,wr=2
195  integer :: aux_stor,band_index,iband,icg,icg_aux,idum
196  integer :: ii,ikg2,ikpt1,ikpt10,ikpt2,ikptsp_prev,inplace,iproc
197  integer :: isppol1,isppol2,istwf10_k,localrdwf
198  integer :: mband1,mband_rd,mband_rw,mcg_aux,me1,me2,mgfft1,mgfft2
199  integer :: my_nspinor1,my_nspinor2
200  integer :: nb_band,nbd1,nbd1_rd,nbd2,nkpt_eff,nproc2,npw1,npw2,nsp
201  integer :: test_cycle,tim_rwwf
202  logical :: out_of_core2
203  character(len=500) :: message
204 !arrays
205  integer,allocatable :: kg1(:,:),kg2_k(:,:),kg_dum(:,:)
206  real(dp) :: kpoint(3),tsec(2)
207  real(dp),allocatable :: cg_aux(:,:),eig_k(:),occ_k(:)
208 
209 ! *************************************************************************
210 
211  call timab(780,1,tsec)
212  call timab(781,1,tsec)
213 
214  icg=0
215 
216 !Init MPI data
217  me1=mpi_enreg1%me_kpt
218  me2=mpi_enreg2%me_kpt
219  nproc2 = mpi_enreg2%nproc_cell
220  out_of_core2=(my_nkpt2/=0.and.mkmem2==0)
221 
222 
223  if((nsppol1==2.and.nspinor2==2).or.(nspinor1==2.and. nsppol2==2))then
224 !  This is not yet possible. See later for a message about where to make the needed modifs.
225 !  EDIT MT 20110707: these modifs are no more needed as they are now done in inwffil
226    write(message, '(a,a,a,a,a,a,a,a,a,i2,a,i2,a,a,i2,a,i2,a,a,a,a)' )ch10,&
227 &   ' newkpt : ERROR -',ch10,&
228 &   '  The wavefunction translator is (still) unable to interchange',ch10,&
229 &   '  spin-polarized wfs and spinor wfs. However,',ch10,&
230 &   '  the input  variables are nsppol1=',nsppol1,', and nspinor1=',nspinor1,ch10,&
231 &   '  the output variables are nsppol2=',nsppol2,', and nspinor2=',nspinor2,ch10,&
232 &   '  Action : use a non-spin-polarized wf to start a spinor wf,',ch10,&
233 &   '           and a non-spinor wf to start a spin-polarized wf.'
234    MSG_ERROR(message)
235  end if
236 
237  my_nspinor1=max(1,nspinor1/mpi_enreg1%nproc_spinor)
238  my_nspinor2=max(1,nspinor2/mpi_enreg2%nproc_spinor)
239  mband1=maxval(nband1(1:nkpt1*nsppol1))
240 
241  if(mkmem1==0 .and. out_of_core2)then
242    mband_rd=min(mband1,(mband2/nspinor2)*nspinor1)
243    if(mcg<mpw1*my_nspinor1*mband_rd)then
244      write(message,'(2(a,i0))')' The dimension mcg= ',mcg,', should be larger than mband_rd= ',mband_rd
245      MSG_BUG(message)
246    end if
247    if(mcg<mband2*mpw2*my_nspinor2)then
248      write(message,'(a,i0,a,a,a,i0,a,i0,a,i2)' )&
249 &     '  The dimension mcg=',mcg,', should be larger than',ch10,&
250 &     '  the product of mband2=',mband2,', mpw2=',mpw2,', and nspinor2=',my_nspinor2
251      MSG_BUG(message)
252    end if
253  end if
254 
255  idum=init_random
256  ikpt10 = 0
257  istwf10_k=0
258  band_index=0
259  icg=0
260 
261  nkpt_eff=nkpt2
262  if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max ) nkpt_eff=nkpt_max
263 
264  mgfft1=maxval(ngfft1(1:3))
265  mgfft2=maxval(ngfft2(1:3))
266  ABI_ALLOCATE(kg1,(3,mpw1))
267  ABI_ALLOCATE(kg2_k,(3,mpw2))
268  ABI_ALLOCATE(kg_dum,(3,0))
269 
270  if (debug>0) then
271    if (me1==0) then
272      write(std_out,'(a)' ) ' newkpt:  kptns1'
273      call prmat (kptns1, 3, nkpt1, 3)
274    end if
275    if (me2==0) then
276      write(std_out,'(a)' ) ' newkpt:  kptns2'
277      call prmat (kptns2, 3, nkpt2, 3)
278    end if
279  end if
280 
281  ikptsp_prev=0
282 
283  call timab(781,2,tsec)
284 
285 !Do outer loop over spins
286  do isppol2=1,nsppol2
287 
288    if (nsppol2==2 .and. me2==0) then
289      write(std_out,'(a,i5)' ) ' newkpt: spin channel isppol2 = ',isppol2
290    end if
291 
292    if (restart==1 .and. out_of_core2) rewind (unkg2)
293    ikg2=0
294 
295 !  Do loop over new k point set
296    do ikpt2=1,nkpt2
297 
298      call timab(782,1,tsec)
299 
300      nbd2=nband2(ikpt2+(isppol2-1)*nkpt2)
301      npw2=npwarr2(ikpt2)
302 
303      if(restart==1)then
304 
305 !      Announce the treatment of k point ikpt
306        if(ikpt2<=nkpt_eff)then
307 !        This message might be overwritten in parallel
308          write(message, '(a,i6,a,i8,a,i4)' )'P newkpt: treating ',nbd2,' bands with npw=',npw2,' for ikpt=',ikpt2
309 !        This message might be overwritten in parallel
310          if(mpi_enreg2%paralbd==1)then
311            do iproc=0,nproc2-1
312              nb_band=0
313              do iband=1,nbd2
314                if(mpi_enreg2%proc_distrb(ikpt2,iband,isppol2) == iproc)nb_band=nb_band+1
315              end do
316              if(nb_band/=0)then
317                write(message, '(a,i6,a,i8,a,i4,a,i4)' ) &
318 &               'P newkpt: treating ',nb_band,' bands with npw=',npw2,' for ikpt=',ikpt2,' by node ',iproc
319              end if
320            end do
321          end if
322          if(mpi_enreg2%paralbd==0) then
323            write(message, '(a,i6,a,i8,a,i4,a,i4)' )&
324 &           'P newkpt: treating ',nbd2,' bands with npw=',npw2,&
325 &           ' for ikpt=',ikpt2,' by node ',mpi_enreg2%proc_distrb(ikpt2,1,isppol2)
326          end if
327          if(prtvol>0)then
328            call wrtout(iout,message,'COLL')
329          end if
330        end if
331 
332 !      Cut the writing if the limit is reached
333        if(ikpt2==nkpt_eff+1)then
334          if(prtvol>0)then
335            call wrtout(iout,' newkpt: prtvol=0 or 1, do not print more k-points.','COLL')
336          end if
337        end if
338 
339 !      End of restart==1
340      end if
341 
342      test_cycle=0
343      if(proc_distrb_cycle(mpi_enreg2%proc_distrb,ikpt2,1,nbd2,isppol2,me2)) test_cycle=1
344      if(test_cycle==1)then
345        if(formeig==0)then
346          eigen(1+band_index : nbd2+band_index) = zero
347 !        occ(1+band_index : nbd2+band_index) = zero
348          band_index=band_index+nbd2
349        else
350          eigen(1+band_index : 2*nbd2**2+band_index) = 0.0_dp
351          band_index=band_index+2*nbd2**2
352        end if
353 !      In the case this k point does not belong to me, cycle
354        if (my_nkpt2==0) cycle
355        if ((mkmem1==0) .and. (ireadwf==1) .and. (mpi_enreg2%paralbd==1))then
356          call WffReadSkipK(formeig,headform1,ikpt2,isppol2,mpi_enreg2,wffinp)
357          ikptsp_prev=ikptsp_prev+1
358        end if
359        cycle
360      end if
361 
362      if(restart==1)then
363 
364        if(mkmem2/=0)then
365          kg2_k(:,1:npw2)=kg2(:,1+ikg2:npw2+ikg2)
366        else if(mkmem2==0)then
367 !        Read the first line of a block and performs some checks on the unkg file.
368          nsp=nspinor2
369          call rdnpw(ikpt2,isppol2,nbd2,npw2,nsp,0,unkg2)
370 !        Read k+g data
371          read (unkg2) kg2_k(1:3,1:npw2)
372        end if
373 
374      end if
375 
376 !    Get ikpt1, the closest k from original set, from indkk
377      ikpt1=indkk(ikpt2,1)
378      if(sppoldbl==2 .and. isppol2==2)ikpt1=indkk(ikpt2+nkpt2,1)
379 
380      npw1=npwarr1(ikpt1)
381      kpoint(:)=kptns1(:,ikpt1)
382 
383 !    Determine the spin polarization of the input data
384      isppol1=isppol2
385      if(nsppol2==2 .and. nsppol1==1)isppol1=1
386 
387      if(restart==2)then
388        if(ikpt2<=nkpt_eff)then
389          write(message,'(a,i4,i8,a,i4,i8)')'- newkpt: read input wf with ikpt,npw=',ikpt1,npw1,', make ikpt,npw=',ikpt2,npw2
390          call wrtout(std_out,message,'PERS')
391          if(iout/=6 .and. me2==0 .and. prtvol>0)then
392            call wrtout(iout,message,'PERS')
393          end if
394        else if(ikpt2==nkpt_eff+1)then
395          write(message,'(a)')'- newkpt : prtvol=0 or 1, do not print more k-points.'
396          call wrtout(std_out,message,'PERS')
397          if(iout/=6 .and. me2==0 .and. prtvol>0)then
398            call wrtout(iout,message,'PERS')
399          end if
400        end if
401      end if
402 
403 !    Set up the number of bands to be read
404      nbd1=nband1(ikpt1+(isppol1-1)*nkpt1)
405      nbd1_rd=min(nbd1,(nbd2/nspinor2)*nspinor1)
406 
407 !    Check that number of bands is not being increased if fill==0 --if so
408 !    print warning and reset new wf file nband2 to only allowed number
409      if ( nbd2/nspinor2 > nbd1/nspinor1 .and. fill==0) then
410        if(ikpt2<=nkpt_eff)then
411          write(message, '(a,i8,a,i8,a,i8)' )' newkpt: nband2=',nbd2,' < nband1=',nbd1,' => reset nband2 to ',nbd1
412          call wrtout(std_out,message,'PERS')
413        end if
414        nbd2=nbd1
415      end if
416 
417 !    Prepare the reading of the wavefunctions: the correct record is selected
418 !    WARNING : works only for GS - for RF the number of record differs
419      if(restart==2 .and. mkmem1==0)then
420 
421        if(debug>0)then
422          write(message, '(a,a,a,a,i5,a,i5,a,a,i5,a,i5)' ) ch10,&
423 &         ' newkpt : about to call randac',ch10,&
424 &         '  for ikpt1=',ikpt1,', ikpt2=',ikpt2,ch10,&
425 &         '  and isppol1=',isppol1,', isppol2=',isppol2
426          call wrtout(std_out,message,'PERS')
427        end if
428 
429        call randac(debug,headform1,ikptsp_prev,ikpt1,isppol1,nband1,nkpt1,nsppol1,wffinp)
430      end if
431 
432 !    Read the data for nbd2 bands at this k point
433 !    Must decide whether an auxiliary storage is needed
434 !    When mkmem1==0 and mkmem2==0 , the cg array should be large enough ...
435 !    When mkmem1==0 and mkmem2/=0 , each k-point block in cg might not be large enough
436 !    however, will read at most (nbd2/nspinor2)*nspinor1 bands from disk
437 !    When mkmem1/=0 , it is supposed that each input k-point block is smaller
438 !    than the corresponding output k-point block, so that the input data
439 !    have been placed already in cg, at the k-point location where they are needed
440      aux_stor=0
441      if(mkmem2/=0 .and. mkmem1==0)then
442        mcg_aux=npw1*my_nspinor1*nbd1
443        if(nbd1_rd<nbd1)mcg_aux=npw1*my_nspinor1*nbd1_rd
444        if( mcg_aux > npw2*my_nspinor2*nbd2 )then
445          aux_stor=1 ; icg_aux=0
446          ABI_ALLOCATE(cg_aux,(2,mcg_aux))
447        end if
448      end if
449 
450      mband_rw=max(nbd1_rd,nbd2)
451      ABI_ALLOCATE(eig_k,(mband_rw*(2*mband_rw)**formeig))
452      if(formeig==0) then
453        ABI_ALLOCATE(occ_k,(mband_rw))
454      else
455        ABI_ALLOCATE(occ_k,(0))
456      end if
457 
458      if(mkmem1/=0 .and. ireadwf==1)then
459 !      Checks that nbd1 and nbd1_rd are equal if eig and occ are input
460        if(nbd1/=nbd1_rd)then
461          write(message,'(a,a,a,i6,a,i6)')&
462 &         '  When mkmem1/=0, one must have nbd1=nbd1_rd, while',ch10,&
463 &         '  nbd1=',nbd1,', and nbd1_rd=',nbd1_rd
464          MSG_BUG(message)
465        end if
466 !      Need to put eigenvalues in eig_k, same for occ
467 !      Note use of band_index, since it is assumed that eigen and occ
468 !      already have spin-k point location structure than output.
469        if(formeig==0)then
470          eig_k(1:nbd1_rd)=eigen(1+band_index : nbd1_rd+band_index)
471 !        occ_k(1:nbd1_rd)=occ(1+band_index : nbd1_rd+band_index)
472        else if(formeig==1)then
473 !        The matrix of eigenvalues has size nbd1 ,  that must be equal
474 !        to nbd1_rd in the case mkmem1/=0)
475          eig_k(1:2*nbd1_rd**2)=eigen(1+band_index : 2*nbd1_rd**2+band_index)
476        end if
477      end if
478 
479      call timab(782,2,tsec)
480 
481 !    Must read the wavefunctions if they are not yet in place
482      if(mkmem1==0 .and. ireadwf==1)then
483 
484        if (debug>0 .and. restart==2) then
485          write(message,'(a,i5,a,a,i5,a,i5,a)' ) &
486 &         ' newkpt : about to call rwwf with ikpt1=',ikpt1,ch10,&
487 &         ' and nband(ikpt1)=',nband1(ikpt1),' nbd2=',nbd2,'.'
488          call wrtout(std_out,message,'PERS')
489        end if
490 
491        if(mpi_enreg1%paralbd==0)tim_rwwf=21
492        if(mpi_enreg1%paralbd==1)tim_rwwf=22
493 
494        if(aux_stor==0)then
495          call rwwf(cg,eig_k,formeig,headform1,icg,ikpt1,isppol1,kg_dum,mband_rw,mcg,mpi_enreg1,&
496 &         nbd1_rd,nbd1,npw1,my_nspinor1,occ_k,1,0,tim_rwwf,wffinp)
497        else
498          icg_aux=0
499          call rwwf(cg_aux,eig_k,formeig,headform1,icg_aux,ikpt1,isppol1,kg_dum,mband_rw,mcg_aux,&
500 &         mpi_enreg1,nbd1_rd,nbd1,npw1,my_nspinor1,occ_k,1,0,tim_rwwf,wffinp)
501        end if
502      end if
503 
504      call timab(783,1,tsec)
505 
506      if(formeig==1 .and. nbd2/=nbd1_rd .and. ireadwf==1)then
507 !      Change the storage of eig_k
508        if(nbd1_rd<nbd2)then
509          do iband=nbd1_rd,1,-1
510 !          The factor of two is for complex eigenvalues
511            do ii=2*nbd2,2*nbd1_rd+1,-1
512              eig_k(ii+(iband-1)*2*nbd2)=huge(0.0_dp)/10.0_dp
513            end do
514            do ii=2*nbd1_rd,1,-1
515              eig_k(ii+(iband-1)*2*nbd2)=eig_k(ii+(iband-1)*2*nbd1_rd)
516            end do
517          end do
518        else if(nbd1_rd>nbd2)then
519          do iband=1,nbd2
520 !          The factor of two is for complex eigenvalues
521            do ii=1,2*nbd2
522              eig_k(ii+(iband-1)*2*nbd2)=eig_k(ii+(iband-1)*2*nbd1_rd)
523            end do
524          end do
525        end if
526      end if
527 
528 !    If change nsppol, must adapt the occupation numbers
529 !    if(nsppol1/=nsppol2)then
530 !    occ_k(1:nbd2)=occ_k(1:nbd2)*nsppol1/dbl(nsppol2)
531 !    then
532 
533 !    In case nsppol1=2 and nspinor2=2, one should read
534 !    the other spin-component, and form a spinor wf here, before calling
535 !    wfconv. One should treat eig_k and occ_k as well.
536 !    A similar operation is to be performed when nspino1=2 and nsppol2=2
537 !    EDIT - MT 20110707: the building of the spinor wf is now done in wfffil
538 !    no need to make it here....
539 
540 !    DEBUG
541 !    write(std_out,*)' newkpt: before wfconv'
542 !    write(std_out,*)' newkpt: mkmem2=',mkmem2
543 !    stop
544 !    ENDDEBUG
545 
546      call timab(783,2,tsec)
547      call timab(784,1,tsec)
548 
549 !    Note the use of mband2, while mband is used inside
550 !    write(std_out,*) 'in newkpt,before wfconv,npw1,npw2',npw1,npw2
551      inplace=1
552      if(aux_stor==0)then
553        call wfconv(ceksp2,cg,cg,debug,ecut1,ecut2,ecut2_eff,&
554 &       eig_k,eig_k,exchn2n3d,formeig,gmet1,gmet2,icg,icg,&
555 &       ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,&
556 &       kg1,kg2_k,kptns1,kptns2,mband_rw,mband_rw,mcg,mcg,&
557 &       mpi_enreg1,mpi_enreg2,mpw1,mpw2,nbd1_rd,nbd2,&
558 &       ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,nsym,&
559 &       occ_k,occ_k,optorth,randalg,restart,rprimd,sppoldbl,symrel,tnons)
560      else
561        call wfconv(ceksp2,cg_aux,cg_aux,debug,ecut1,ecut2,ecut2_eff,&
562 &       eig_k,eig_k,exchn2n3d,formeig,gmet1,gmet2,icg_aux,icg_aux,&
563 &       ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,&
564 &       kg1,kg2_k,kptns1,kptns2,mband_rw,mband_rw,mcg,mcg,&
565 &       mpi_enreg1,mpi_enreg2,mpw1,mpw2,nbd1_rd,nbd2,&
566 &       ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,nsym,&
567 &       occ_k,occ_k,optorth,randalg,restart,rprimd,sppoldbl,symrel,tnons)
568      end if
569 
570      call timab(784,2,tsec)
571 
572 !    Finally write new wf to disk file or save in permanent file
573      if(mkmem2==0)then
574 
575 !      Note that in this case, we are sure aux_stor==0
576        if(mpi_enreg2%paralbd==0)tim_rwwf=21
577        if(mpi_enreg2%paralbd==1)tim_rwwf=22
578        call rwwf(cg,eig_k,formeig,0,0,ikpt2,isppol2,kg2_k,nbd2,mcg,mpi_enreg2,&
579 &       nbd2,nbd2,npw2,my_nspinor2,occ_k,wr,1,tim_rwwf,wffout)
580 
581      end if
582 
583      call timab(785,1,tsec)
584 
585      if(mkmem2/=0)then
586        if(aux_stor==1)then
587          cg(:,1+icg:npw2*nbd2*my_nspinor2+icg)=cg_aux(:,1:npw2*nbd2*my_nspinor2)
588          ABI_DEALLOCATE(cg_aux)
589        end if
590 
591        icg=icg+npw2*nbd2*my_nspinor2
592        ikg2=ikg2+npw2
593      end if
594 
595      eigen(1+band_index:nbd2*(2*nbd2)**formeig+band_index) = eig_k(1:nbd2*(2*nbd2)**formeig)
596 !    occ(1+band_index:nbd2+band_index)=occ_k(1:nbd2)
597 
598      if(formeig==0)then
599        band_index=band_index+nbd2
600      else if(formeig==1)then
601        band_index=band_index+2*nbd2**2
602      end if
603 
604      ABI_DEALLOCATE(eig_k)
605      ABI_DEALLOCATE(occ_k)
606 
607      call timab(785,2,tsec)
608 
609    end do ! ikpt2
610  end do ! isppol2
611 
612  call timab(786,1,tsec)
613 
614  if(xmpi_paral==1)then
615 !  Transmit eigenvalues (not yet occupation numbers)
616 !  newkpt.F90 is not yet suited for RF format
617 !  This routine works in both localrdwf=0 or 1 cases.
618 !  However, in the present routine, localrdwf is to be considered
619 !  as 1 always, since the transfer has been made in wfsinp .
620    localrdwf=1
621    call pareigocc(eigen,formeig,localrdwf,mpi_enreg2,mband2,nband2,nkpt2,nsppol2,occ,1)
622  end if
623 
624  ABI_DEALLOCATE(kg1)
625  ABI_DEALLOCATE(kg2_k)
626  ABI_DEALLOCATE(kg_dum)
627 
628  call timab(786,2,tsec)
629  call timab(780,2,tsec)
630 
631 end subroutine newkpt