TABLE OF CONTENTS


ABINIT/m_inwffil [ Modules ]

[ Top ] [ Modules ]

NAME

  m_inwffil

FUNCTION

  Do initialization of wavefunction files.

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, AR, MB, MVer, ZL, MB, TD)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_inwffil
28 
29  use defs_basis
30  use defs_abitypes
31  use defs_wvltypes
32  use m_abicore
33  use m_wffile
34  use m_wfk
35  use m_errors
36  use m_xmpi
37  use m_nctk
38  use m_hdr
39 #if defined HAVE_MPI2
40  use mpi
41 #endif
42 
43  use m_time,     only : timab
44  use m_io_tools, only : file_exists, get_unit
45  use m_geometry, only : getspinrot
46  use m_pptools,  only : prmat
47  use m_symtk,    only : matr3inv, mati3inv
48  use m_cgtools,  only : cg_envlop, pw_orthon
49  use m_fftcore,  only : kpgsph, sphere, sphereboundary
50  use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_io
51  use m_mpinfo,   only : destroy_mpi_enreg, copy_mpi_enreg, proc_distrb_cycle
52  use m_kg,       only : kpgio, ph1d3d, getph
53  use m_kpts,     only : listkk
54  use m_occ,      only : pareigocc
55  use m_rwwf,     only : rwwf, WffReadSkipK
56  use m_wvl_wfsinp, only : wvl_wfsinp_disk, wvl_wfsinp_scratch
57 
58  implicit none
59 
60 #if defined HAVE_MPI1
61  include 'mpif.h'
62 #endif
63 
64  private

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.

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=information about MPI parallelization, for the input wf file
  mpi_enreg2=information 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

2104 subroutine newkpt(ceksp2,cg,debug,ecut1,ecut2,ecut2_eff,eigen,exchn2n3d,fill,&
2105 &                  formeig,gmet1,gmet2,headform1,indkk,iout,ireadwf,&
2106 &                  istwfk1,istwfk2,kg2,kptns1,kptns2,mband2,mcg,mkmem1,mkmem2,&
2107 &                  mpi_enreg1,mpi_enreg2,mpw1,mpw2,my_nkpt2,nband1,nband2,&
2108 &                  ngfft1,ngfft2,nkpt1,nkpt2,npwarr1,npwarr2,nspinor1,nspinor2,&
2109 &                  nsppol1,nsppol2,nsym,occ,optorth,prtvol,randalg,restart,rprimd,&
2110 &                  sppoldbl,symrel,tnons,unkg2,wffinp,wffout)
2111 
2112 
2113 !This section has been created automatically by the script Abilint (TD).
2114 !Do not modify the following lines by hand.
2115 #undef ABI_FUNC
2116 #define ABI_FUNC 'newkpt'
2117 !End of the abilint section
2118 
2119  implicit none
2120 
2121 !Arguments ------------------------------------
2122 !scalars
2123  integer,intent(in) :: ceksp2,debug,exchn2n3d,fill,formeig,headform1,iout
2124  integer,intent(in) :: ireadwf,mband2,mcg,mkmem1,mkmem2,mpw1,mpw2,my_nkpt2,nkpt1,nkpt2
2125  integer,intent(in) :: nspinor1,nspinor2,nsppol1,nsppol2,nsym,optorth,prtvol,restart
2126  integer,intent(in) :: randalg,sppoldbl,unkg2
2127  real(dp),intent(in) :: ecut1,ecut2,ecut2_eff
2128  type(MPI_type),intent(inout) :: mpi_enreg1,mpi_enreg2
2129  type(wffile_type),intent(inout) :: wffinp,wffout
2130 !arrays
2131  integer,intent(in) :: indkk(nkpt2*sppoldbl,6),istwfk1(nkpt1),istwfk2(nkpt2)
2132  integer,intent(in) :: kg2(3,mpw2*mkmem2),nband1(nkpt1*nsppol1)
2133  integer,intent(in) :: nband2(nkpt2*nsppol2),ngfft1(18),ngfft2(18)
2134  integer,intent(in) :: npwarr1(nkpt1),npwarr2(nkpt2),symrel(3,3,nsym)
2135  real(dp),intent(in) :: gmet1(3,3),gmet2(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2)
2136  real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym)
2137  real(dp),intent(inout) :: cg(2,mcg) !vz_i pw_orthon vecnm
2138  real(dp),intent(inout) :: eigen(mband2*(2*mband2)**formeig*nkpt2*nsppol2)!vz_i newocc
2139  real(dp),intent(inout) :: occ(mband2*nkpt2*nsppol2) !vz_i
2140 
2141 !Local variables-------------------------------
2142 !scalars
2143  integer,parameter :: init_random=-5,nkpt_max=50,tobox=1,tosph=-1,wr=2
2144  integer :: aux_stor,band_index,iband,icg,icg_aux,idum
2145  integer :: ii,ikg2,ikpt1,ikpt10,ikpt2,ikptsp_prev,inplace,iproc
2146  integer :: isppol1,isppol2,istwf10_k,localrdwf
2147  integer :: mband1,mband_rd,mband_rw,mcg_aux,me1,me2,mgfft1,mgfft2
2148  integer :: my_nspinor1,my_nspinor2
2149  integer :: nb_band,nbd1,nbd1_rd,nbd2,nkpt_eff,nproc2,npw1,npw2,nsp
2150  integer :: test_cycle,tim_rwwf
2151  logical :: out_of_core2
2152  character(len=500) :: message
2153 !arrays
2154  integer,allocatable :: kg1(:,:),kg2_k(:,:),kg_dum(:,:)
2155  real(dp) :: kpoint(3),tsec(2)
2156  real(dp),allocatable :: cg_aux(:,:),eig_k(:),occ_k(:)
2157 
2158 ! *************************************************************************
2159 
2160  call timab(780,1,tsec)
2161  call timab(781,1,tsec)
2162 
2163  icg=0
2164 
2165 !Init MPI data
2166  me1=mpi_enreg1%me_kpt
2167  me2=mpi_enreg2%me_kpt
2168  nproc2 = mpi_enreg2%nproc_cell
2169  out_of_core2=(my_nkpt2/=0.and.mkmem2==0)
2170 
2171 
2172  if((nsppol1==2.and.nspinor2==2).or.(nspinor1==2.and. nsppol2==2))then
2173 !  This is not yet possible. See later for a message about where to make the needed modifs.
2174 !  EDIT MT 20110707: these modifs are no more needed as they are now done in inwffil
2175    write(message, '(5a,i2,a,i2,2a,i2,a,i2,4a)' ) &
2176 &   'The wavefunction translator is (still) unable to interchange',ch10,&
2177 &   'spin-polarized wfs and spinor wfs. However,',ch10,&
2178 &   'the input  variables are nsppol1=',nsppol1,', and nspinor1=',nspinor1,ch10,&
2179 &   'the output variables are nsppol2=',nsppol2,', and nspinor2=',nspinor2,ch10,&
2180 &   'Action: use a non-spin-polarized wf to start a spinor wf,',ch10,&
2181 &   '        and a non-spinor wf to start a spin-polarized wf.'
2182    MSG_ERROR(message)
2183  end if
2184 
2185  my_nspinor1=max(1,nspinor1/mpi_enreg1%nproc_spinor)
2186  my_nspinor2=max(1,nspinor2/mpi_enreg2%nproc_spinor)
2187  mband1=maxval(nband1(1:nkpt1*nsppol1))
2188 
2189  if(mkmem1==0 .and. out_of_core2)then
2190    mband_rd=min(mband1,(mband2/nspinor2)*nspinor1)
2191    if(mcg<mpw1*my_nspinor1*mband_rd)then
2192      write(message,'(2(a,i0))')' The dimension mcg= ',mcg,', should be larger than mband_rd= ',mband_rd
2193      MSG_BUG(message)
2194    end if
2195    if(mcg<mband2*mpw2*my_nspinor2)then
2196      write(message,'(a,i0,a,a,a,i0,a,i0,a,i2)' )&
2197 &     'The dimension mcg= ',mcg,', should be larger than',ch10,&
2198 &     'the product of mband2= ',mband2,', mpw2= ',mpw2,', and nspinor2= ',my_nspinor2
2199      MSG_BUG(message)
2200    end if
2201  end if
2202 
2203  idum=init_random
2204  ikpt10 = 0
2205  istwf10_k=0
2206  band_index=0
2207  icg=0
2208 
2209  nkpt_eff=nkpt2
2210  if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max ) nkpt_eff=nkpt_max
2211 
2212  mgfft1=maxval(ngfft1(1:3))
2213  mgfft2=maxval(ngfft2(1:3))
2214  ABI_ALLOCATE(kg1,(3,mpw1))
2215  ABI_ALLOCATE(kg2_k,(3,mpw2))
2216  ABI_ALLOCATE(kg_dum,(3,0))
2217 
2218  if (debug>0) then
2219    if (me1==0) then
2220      write(std_out,'(a)' ) ' newkpt:  kptns1'
2221      call prmat (kptns1, 3, nkpt1, 3)
2222    end if
2223    if (me2==0) then
2224      write(std_out,'(a)' ) ' newkpt:  kptns2'
2225      call prmat (kptns2, 3, nkpt2, 3)
2226    end if
2227  end if
2228 
2229  ikptsp_prev=0
2230 
2231  call timab(781,2,tsec)
2232 
2233 !Do outer loop over spins
2234  do isppol2=1,nsppol2
2235 
2236    if (nsppol2==2 .and. me2==0) then
2237      write(std_out,'(a,i5)' ) ' newkpt: spin channel isppol2 = ',isppol2
2238    end if
2239 
2240    if (restart==1 .and. out_of_core2) rewind (unkg2)
2241    ikg2=0
2242 
2243 !  Do loop over new k point set
2244    do ikpt2=1,nkpt2
2245 
2246      call timab(782,1,tsec)
2247 
2248      nbd2=nband2(ikpt2+(isppol2-1)*nkpt2)
2249      npw2=npwarr2(ikpt2)
2250 
2251      if(restart==1)then
2252 
2253 !      Announce the treatment of k point ikpt
2254        if(ikpt2<=nkpt_eff)then
2255 !        This message might be overwritten in parallel
2256          write(message, '(a,i6,a,i8,a,i4)' )'P newkpt: treating ',nbd2,' bands with npw=',npw2,' for ikpt=',ikpt2
2257 !        This message might be overwritten in parallel
2258          if(mpi_enreg2%paralbd==1)then
2259            do iproc=0,nproc2-1
2260              nb_band=0
2261              do iband=1,nbd2
2262                if(mpi_enreg2%proc_distrb(ikpt2,iband,isppol2) == iproc)nb_band=nb_band+1
2263              end do
2264              if(nb_band/=0)then
2265                write(message, '(a,i6,a,i8,a,i4,a,i4)' ) &
2266 &               'P newkpt: treating ',nb_band,' bands with npw=',npw2,' for ikpt=',ikpt2,' by node ',iproc
2267              end if
2268            end do
2269          end if
2270          if(mpi_enreg2%paralbd==0) then
2271            write(message, '(a,i6,a,i8,a,i4,a,i4)' )&
2272 &           'P newkpt: treating ',nbd2,' bands with npw=',npw2,&
2273 &           ' for ikpt=',ikpt2,' by node ',mpi_enreg2%proc_distrb(ikpt2,1,isppol2)
2274          end if
2275          if(prtvol>0)then
2276            call wrtout(iout,message,'COLL')
2277          end if
2278        end if
2279 
2280 !      Cut the writing if the limit is reached
2281        if(ikpt2==nkpt_eff+1)then
2282          if(prtvol>0)then
2283            call wrtout(iout,' newkpt: prtvol=0 or 1, do not print more k-points.','COLL')
2284          end if
2285        end if
2286 
2287 !      End of restart==1
2288      end if
2289 
2290      test_cycle=0
2291      if(proc_distrb_cycle(mpi_enreg2%proc_distrb,ikpt2,1,nbd2,isppol2,me2)) test_cycle=1
2292      if(test_cycle==1)then
2293        if(formeig==0)then
2294          eigen(1+band_index : nbd2+band_index) = zero
2295 !        occ(1+band_index : nbd2+band_index) = zero
2296          band_index=band_index+nbd2
2297        else
2298          eigen(1+band_index : 2*nbd2**2+band_index) = 0.0_dp
2299          band_index=band_index+2*nbd2**2
2300        end if
2301 !      In the case this k point does not belong to me, cycle
2302        if (my_nkpt2==0) cycle
2303        if ((mkmem1==0) .and. (ireadwf==1) .and. (mpi_enreg2%paralbd==1))then
2304          call WffReadSkipK(formeig,headform1,ikpt2,isppol2,mpi_enreg2,wffinp)
2305          ikptsp_prev=ikptsp_prev+1
2306        end if
2307        cycle
2308      end if
2309 
2310      if(restart==1)then
2311 
2312        if(mkmem2/=0)then
2313          kg2_k(:,1:npw2)=kg2(:,1+ikg2:npw2+ikg2)
2314        else if(mkmem2==0)then
2315 !        Read the first line of a block and performs some checks on the unkg file.
2316          MSG_ERROR("mkmem2 == 0 and rdnpw are not supported anymore.")
2317          nsp=nspinor2
2318          !call rdnpw(ikpt2,isppol2,nbd2,npw2,nsp,0,unkg2)
2319 !        Read k+g data
2320          read (unkg2) kg2_k(1:3,1:npw2)
2321        end if
2322 
2323      end if
2324 
2325 !    Get ikpt1, the closest k from original set, from indkk
2326      ikpt1=indkk(ikpt2,1)
2327      if(sppoldbl==2 .and. isppol2==2)ikpt1=indkk(ikpt2+nkpt2,1)
2328 
2329      npw1=npwarr1(ikpt1)
2330      kpoint(:)=kptns1(:,ikpt1)
2331 
2332 !    Determine the spin polarization of the input data
2333      isppol1=isppol2
2334      if(nsppol2==2 .and. nsppol1==1)isppol1=1
2335 
2336      if(restart==2)then
2337        if(ikpt2<=nkpt_eff)then
2338          write(message,'(a,i4,i8,a,i4,i8)')'- newkpt: read input wf with ikpt,npw=',ikpt1,npw1,', make ikpt,npw=',ikpt2,npw2
2339          call wrtout(std_out,message,'PERS')
2340          if(iout/=6 .and. me2==0 .and. prtvol>0)then
2341            call wrtout(iout,message,'PERS')
2342          end if
2343        else if(ikpt2==nkpt_eff+1)then
2344          call wrtout(std_out, '- newkpt : prtvol=0 or 1, do not print more k-points.', 'PERS')
2345          if(iout/=6 .and. me2==0 .and. prtvol>0)then
2346            call wrtout(iout,message,'PERS')
2347          end if
2348        end if
2349      end if
2350 
2351 !    Set up the number of bands to be read
2352      nbd1=nband1(ikpt1+(isppol1-1)*nkpt1)
2353      nbd1_rd=min(nbd1,(nbd2/nspinor2)*nspinor1)
2354 
2355 !    Check that number of bands is not being increased if fill==0 --if so
2356 !    print warning and reset new wf file nband2 to only allowed number
2357      if ( nbd2/nspinor2 > nbd1/nspinor1 .and. fill==0) then
2358        if(ikpt2<=nkpt_eff)then
2359          write(message, '(a,i8,a,i8,a,i8)' )' newkpt: nband2=',nbd2,' < nband1=',nbd1,' => reset nband2 to ',nbd1
2360          call wrtout(std_out,message,'PERS')
2361        end if
2362        nbd2=nbd1
2363      end if
2364 
2365 !    Prepare the reading of the wavefunctions: the correct record is selected
2366 !    WARNING : works only for GS - for RF the number of record differs
2367      if(restart==2 .and. mkmem1==0)then
2368        MSG_ERROR("mkmem1 == 0 has been removed.")
2369 
2370        if(debug>0)then
2371          write(message, '(a,a,a,a,i5,a,i5,a,a,i5,a,i5)' ) ch10,&
2372          ' newkpt : about to call randac',ch10,&
2373          '  for ikpt1=',ikpt1,', ikpt2=',ikpt2,ch10,&
2374          '  and isppol1=',isppol1,', isppol2=',isppol2
2375          call wrtout(std_out,message,'PERS')
2376        end if
2377 
2378        !call randac(debug,headform1,ikptsp_prev,ikpt1,isppol1,nband1,nkpt1,nsppol1,wffinp)
2379      end if
2380 
2381 !    Read the data for nbd2 bands at this k point
2382 !    Must decide whether an auxiliary storage is needed
2383 !    When mkmem1==0 and mkmem2==0 , the cg array should be large enough ...
2384 !    When mkmem1==0 and mkmem2/=0 , each k-point block in cg might not be large enough
2385 !    however, will read at most (nbd2/nspinor2)*nspinor1 bands from disk
2386 !    When mkmem1/=0 , it is supposed that each input k-point block is smaller
2387 !    than the corresponding output k-point block, so that the input data
2388 !    have been placed already in cg, at the k-point location where they are needed
2389      aux_stor=0
2390      if(mkmem2/=0 .and. mkmem1==0)then
2391        mcg_aux=npw1*my_nspinor1*nbd1
2392        if(nbd1_rd<nbd1)mcg_aux=npw1*my_nspinor1*nbd1_rd
2393        if( mcg_aux > npw2*my_nspinor2*nbd2 )then
2394          aux_stor=1 ; icg_aux=0
2395          ABI_ALLOCATE(cg_aux,(2,mcg_aux))
2396        end if
2397      end if
2398 
2399      mband_rw=max(nbd1_rd,nbd2)
2400      ABI_ALLOCATE(eig_k,(mband_rw*(2*mband_rw)**formeig))
2401      if(formeig==0) then
2402        ABI_ALLOCATE(occ_k,(mband_rw))
2403      else
2404        ABI_ALLOCATE(occ_k,(0))
2405      end if
2406 
2407      if(mkmem1/=0 .and. ireadwf==1)then
2408 !      Checks that nbd1 and nbd1_rd are equal if eig and occ are input
2409        if(nbd1/=nbd1_rd)then
2410          write(message,'(a,a,a,i6,a,i6)')&
2411 &         'When mkmem1/=0, one must have nbd1=nbd1_rd, while',ch10,&
2412 &         'nbd1 = ',nbd1,', and nbd1_rd = ',nbd1_rd
2413          MSG_BUG(message)
2414        end if
2415 !      Need to put eigenvalues in eig_k, same for occ
2416 !      Note use of band_index, since it is assumed that eigen and occ
2417 !      already have spin-k point location structure than output.
2418        if(formeig==0)then
2419          eig_k(1:nbd1_rd)=eigen(1+band_index : nbd1_rd+band_index)
2420 !        occ_k(1:nbd1_rd)=occ(1+band_index : nbd1_rd+band_index)
2421        else if(formeig==1)then
2422 !        The matrix of eigenvalues has size nbd1 ,  that must be equal
2423 !        to nbd1_rd in the case mkmem1/=0)
2424          eig_k(1:2*nbd1_rd**2)=eigen(1+band_index : 2*nbd1_rd**2+band_index)
2425        end if
2426      end if
2427 
2428      call timab(782,2,tsec)
2429 
2430 !    Must read the wavefunctions if they are not yet in place
2431      if(mkmem1==0 .and. ireadwf==1)then
2432 
2433        if (debug>0 .and. restart==2) then
2434          write(message,'(a,i5,a,a,i5,a,i5,a)' ) &
2435 &         ' newkpt : about to call rwwf with ikpt1=',ikpt1,ch10,&
2436 &         ' and nband(ikpt1)=',nband1(ikpt1),' nbd2=',nbd2,'.'
2437          call wrtout(std_out,message,'PERS')
2438        end if
2439 
2440        if(mpi_enreg1%paralbd==0)tim_rwwf=21
2441        if(mpi_enreg1%paralbd==1)tim_rwwf=22
2442 
2443        if(aux_stor==0)then
2444          call rwwf(cg,eig_k,formeig,headform1,icg,ikpt1,isppol1,kg_dum,mband_rw,mcg,mpi_enreg1,&
2445 &         nbd1_rd,nbd1,npw1,my_nspinor1,occ_k,1,0,tim_rwwf,wffinp)
2446        else
2447          icg_aux=0
2448          call rwwf(cg_aux,eig_k,formeig,headform1,icg_aux,ikpt1,isppol1,kg_dum,mband_rw,mcg_aux,&
2449 &         mpi_enreg1,nbd1_rd,nbd1,npw1,my_nspinor1,occ_k,1,0,tim_rwwf,wffinp)
2450        end if
2451      end if
2452 
2453      call timab(783,1,tsec)
2454 
2455      if(formeig==1 .and. nbd2/=nbd1_rd .and. ireadwf==1)then
2456 !      Change the storage of eig_k
2457        if(nbd1_rd<nbd2)then
2458          do iband=nbd1_rd,1,-1
2459 !          The factor of two is for complex eigenvalues
2460            do ii=2*nbd2,2*nbd1_rd+1,-1
2461              eig_k(ii+(iband-1)*2*nbd2)=huge(0.0_dp)/10.0_dp
2462            end do
2463            do ii=2*nbd1_rd,1,-1
2464              eig_k(ii+(iband-1)*2*nbd2)=eig_k(ii+(iband-1)*2*nbd1_rd)
2465            end do
2466          end do
2467        else if(nbd1_rd>nbd2)then
2468          do iband=1,nbd2
2469 !          The factor of two is for complex eigenvalues
2470            do ii=1,2*nbd2
2471              eig_k(ii+(iband-1)*2*nbd2)=eig_k(ii+(iband-1)*2*nbd1_rd)
2472            end do
2473          end do
2474        end if
2475      end if
2476 
2477 !    If change nsppol, must adapt the occupation numbers
2478 !    if(nsppol1/=nsppol2)then
2479 !    occ_k(1:nbd2)=occ_k(1:nbd2)*nsppol1/dbl(nsppol2)
2480 !    then
2481 
2482 !    In case nsppol1=2 and nspinor2=2, one should read
2483 !    the other spin-component, and form a spinor wf here, before calling
2484 !    wfconv. One should treat eig_k and occ_k as well.
2485 !    A similar operation is to be performed when nspino1=2 and nsppol2=2
2486 !    EDIT - MT 20110707: the building of the spinor wf is now done in wfffil
2487 !    no need to make it here....
2488 
2489 !    DEBUG
2490 !    write(std_out,*)' newkpt: before wfconv'
2491 !    write(std_out,*)' newkpt: mkmem2=',mkmem2
2492 !    stop
2493 !    ENDDEBUG
2494 
2495      call timab(783,2,tsec)
2496      call timab(784,1,tsec)
2497 
2498 !    Note the use of mband2, while mband is used inside
2499 !    write(std_out,*) 'in newkpt,before wfconv,npw1,npw2',npw1,npw2
2500      inplace=1
2501      if(aux_stor==0)then
2502        call wfconv(ceksp2,cg,cg,debug,ecut1,ecut2,ecut2_eff,&
2503 &       eig_k,eig_k,exchn2n3d,formeig,gmet1,gmet2,icg,icg,&
2504 &       ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,&
2505 &       kg1,kg2_k,kptns1,kptns2,mband_rw,mband_rw,mcg,mcg,&
2506 &       mpi_enreg1,mpi_enreg2,mpw1,mpw2,nbd1_rd,nbd2,&
2507 &       ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,nsym,&
2508 &       occ_k,occ_k,optorth,randalg,restart,rprimd,sppoldbl,symrel,tnons)
2509      else
2510        call wfconv(ceksp2,cg_aux,cg_aux,debug,ecut1,ecut2,ecut2_eff,&
2511 &       eig_k,eig_k,exchn2n3d,formeig,gmet1,gmet2,icg_aux,icg_aux,&
2512 &       ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,&
2513 &       kg1,kg2_k,kptns1,kptns2,mband_rw,mband_rw,mcg,mcg,&
2514 &       mpi_enreg1,mpi_enreg2,mpw1,mpw2,nbd1_rd,nbd2,&
2515 &       ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,nsym,&
2516 &       occ_k,occ_k,optorth,randalg,restart,rprimd,sppoldbl,symrel,tnons)
2517      end if
2518 
2519      call timab(784,2,tsec)
2520 
2521 !    Finally write new wf to disk file or save in permanent file
2522      if(mkmem2==0)then
2523 
2524 !      Note that in this case, we are sure aux_stor==0
2525        if(mpi_enreg2%paralbd==0)tim_rwwf=21
2526        if(mpi_enreg2%paralbd==1)tim_rwwf=22
2527        call rwwf(cg,eig_k,formeig,0,0,ikpt2,isppol2,kg2_k,nbd2,mcg,mpi_enreg2,&
2528 &       nbd2,nbd2,npw2,my_nspinor2,occ_k,wr,1,tim_rwwf,wffout)
2529 
2530      end if
2531 
2532      call timab(785,1,tsec)
2533 
2534      if(mkmem2/=0)then
2535        if(aux_stor==1)then
2536          cg(:,1+icg:npw2*nbd2*my_nspinor2+icg)=cg_aux(:,1:npw2*nbd2*my_nspinor2)
2537          ABI_DEALLOCATE(cg_aux)
2538        end if
2539 
2540        icg=icg+npw2*nbd2*my_nspinor2
2541        ikg2=ikg2+npw2
2542      end if
2543 
2544      eigen(1+band_index:nbd2*(2*nbd2)**formeig+band_index) = eig_k(1:nbd2*(2*nbd2)**formeig)
2545 !    occ(1+band_index:nbd2+band_index)=occ_k(1:nbd2)
2546 
2547      if(formeig==0)then
2548        band_index=band_index+nbd2
2549      else if(formeig==1)then
2550        band_index=band_index+2*nbd2**2
2551      end if
2552 
2553      ABI_DEALLOCATE(eig_k)
2554      ABI_DEALLOCATE(occ_k)
2555 
2556      call timab(785,2,tsec)
2557 
2558    end do ! ikpt2
2559  end do ! isppol2
2560 
2561  call timab(786,1,tsec)
2562 
2563  if(xmpi_paral==1)then
2564 !  Transmit eigenvalues (not yet occupation numbers)
2565 !  newkpt.F90 is not yet suited for RF format
2566 !  This routine works in both localrdwf=0 or 1 cases.
2567 !  However, in the present routine, localrdwf is to be considered
2568 !  as 1 always, since the transfer has been made in wfsinp .
2569    localrdwf=1
2570    call pareigocc(eigen,formeig,localrdwf,mpi_enreg2,mband2,nband2,nkpt2,nsppol2,occ,1)
2571  end if
2572 
2573  ABI_DEALLOCATE(kg1)
2574  ABI_DEALLOCATE(kg2_k)
2575  ABI_DEALLOCATE(kg_dum)
2576 
2577  call timab(786,2,tsec)
2578  call timab(780,2,tsec)
2579 
2580 end subroutine newkpt

ABINIT/wfconv [ Functions ]

[ Top ] [ Functions ]

NAME

 wfconv

FUNCTION

 This subroutine treats the wavefunctions for one k point,
 and converts them to other parameters.

INPUTS

  ceksp2=if 1, center the output sphere of pw on Gamma; if 0, on each k-point (usual).
  cg1(2,mcg1)=wavefunction array
  debug= if 1, print some messages ; otherwise, 0.
  ecut1=kinetic energy cutoffs for basis sphere 1 (hartree)
  ecut2=kinetic energy cutoff beyond which the coefficients of wf2 vanish (Ha)
  ecut2_eff=kinetic energy cut-off for basis sphere 2 (hartree)
  eig_k1(mband1*(2*mband1)**formeig)=eigenvalues
  exchn2n3d=if 1, n2 and n3 are exchanged
  formeig option (format of the eigenvalues and eigenvector) :
   0 => ground-state format (initialisation of
        eigenvectors with random numbers, vector of eigenvalues)
   1 => respfn format (initialisation of
        eigenvectors with 0 s, hermitian matrix of eigenvalues)
  gmet1(3,3)=reciprocal space metric (bohr^-2) for input wf
  gmet2(3,3)=reciprocal space metric (bohr^-2) for output wf
  icg1=shift to be given to the location of the data in the array cg1
  icg2=shift to be given to the location of the data in the array cg2
  ikpt1=number of the k point actually treated (input wf numbering)
  ikpt10=number of the k point previously treated (input wf numbering)
  ikpt2=number of the k point actually treated (output numbering)
  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 kpt2.
   indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
  inplace= if 0, cg1 and cg2 are different in the calling routine,
           if 1, cg1 and cg2 are identical (they have the same memory location)
    This is also true for the pairs (eig_k1,eig_k2) and (occ_k1,occ_k2)
  isppol2=spin variable for output wavefunctions
  istwfk1(nkpt1)=input parameter that describes the storage of wfs in set1
  istwfk2(nkpt2)=input parameter that describes the storage of wfs in set2
  kg1(3,mpw1)=dimensionless coords of G vecs in basis sphere at k point (input wf)
  kg2(3,mpw2)=dimensionless coords of G vecs in basis sphere at k point (output wf)
  kptns1(3,nkpt1)=k point set for input wavefunctions
  kptns2(3,nkpt2)=k point set for output wavefunctions
  mband1=dimension of eig_k1 and occ_k1 arrays
  mband2=dimension of eig_k2 and occ_k2 arrays
  mcg1=dimension of cg1 array (at least npw1*nspinor1*nbd1)
  mcg2=dimension of cg2 array (at least npw2*nspinor2*nbd2)
  mpi_enreg1=information about MPI parallelization for set 1
  mpi_enreg2=information about MPI parallelization for set 2
  mpw1=dimension of kg1, can be set to 0 if not needed
  mpw2=dimension of kg2, can be set to 0 if not needed
  nbd1=number of bands contained in cg1,eig_k1,occ_k1 at this k-point - spin (at input)
  nbd2=number of bands contained in cg2,eig_k2,occ_k2 at this k-point - spin (at output)
  ngfft1(18)=all needed information about 3D FFT, for input wavefunctions
  ngfft2(18)=all needed information about 3D FFT, for output wavefunctions
             see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt1=number of k points for input wavefunctions
  nkpt2=number of k points for output wavefunctions
  npw1=number of planewaves for input wavefunctions
  npw2=number of planewaves for output wavefunctions
  nspinor1=number of spinors for input wavefunctions
  nspinor2=number of spinors for output wavefunctions
  nsym=number of symmetry elements in space group
  occ_k1(mband1)=occupation numbers
  optorth=1 if the WFs are orthogonalized before leaving the routine
  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)
  rprimd2(3,3)=dimensional primitive translations for real space (bohr)
   needed only for the spinor rotation
  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

OUTPUT

  cg2(2,mcg2)=wavefunction array
  eig_k2(mband2*(2*mband2)**formeig)=eigenvalues
  occ_k2(mband2)=occupation (completed with zeros)

SIDE EFFECTS

 Input/Output:
  ikpt10=at input, number of the k point previously treated (input wf numbering)
     (if this is the first call for the present k point set, ikpt10 should be 0)
         at output, number of the k point just treated (input wf numbering)
  kg1, kg2, npw1 and npw2 should not be modified by kpgsph (TD).

NOTES

 Note that this routine can make an in-place conversion
 (see the input variable "inplace"),
 if cg1 and cg2 are equal, as well as the pairs (icg1,icg2),
 (eig_k1,eig_k2),(occ_k1,occ_k2) and (mband1,mband2)

 It can also be used to fill or to initialize wavefunctions
 at one k point
 (filling with random numbers or 0''s, according to the value
 of formeig), if the input number of bands (nbd1) is 0.
 In the latter case, one should use the same values of input
 wavefunction parameters
 than for output wavefunction parameters, except nbd1.

 The input parameters are indexed with 1, the output parameters
 are indexed with 2.

 Some of the arguments are arrays dimensioned with nkpt1 or nkpt2.
 Note that for these, only the elements for ikpt1 or ikpt2 will be used.

 The number of input bands must already be minimal at the input.
 This means, when input and output nspinor are equal : nbd1<nbd2
 When the two nspinor differ, one must have nbd1/nspinor1<nbd2/nspinor2

PARENTS

      newkpt,wfsinp

CHILDREN

      cg_envlop,getph,getspinrot,kpgsph,mati3inv,ph1d3d,pw_orthon,sphere
      sphereboundary,timab,wrtout,xmpi_sum

SOURCE

2708 subroutine wfconv(ceksp2,cg1,cg2,debug,ecut1,ecut2,ecut2_eff,&
2709 & eig_k1,eig_k2,exchn2n3d,formeig,gmet1,gmet2,icg1,icg2,&
2710 & ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,&
2711 & kg1,kg2,kptns1,kptns2,mband1,mband2,mcg1,mcg2,mpi_enreg1,mpi_enreg2,&
2712 & mpw1,mpw2,nbd1,nbd2,ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,&
2713 & nsym,occ_k1,occ_k2,optorth,randalg,restart,rprimd2,sppoldbl,symrel,tnons)
2714 
2715 
2716 !This section has been created automatically by the script Abilint (TD).
2717 !Do not modify the following lines by hand.
2718 #undef ABI_FUNC
2719 #define ABI_FUNC 'wfconv'
2720 !End of the abilint section
2721 
2722  implicit none
2723 
2724 !Arguments ------------------------------------
2725 !scalars
2726  integer,intent(in) :: ceksp2,debug,exchn2n3d,formeig,icg1,icg2,ikpt1
2727  integer,intent(in) :: ikpt2,inplace,isppol2,mband1,mband2,mcg1,mcg2,mpw1,mpw2
2728  integer,intent(in) :: nbd1,nbd2,nkpt1,nkpt2,nspinor1,nspinor2,nsym
2729  integer,intent(in) :: optorth,randalg,restart,sppoldbl
2730  integer,intent(inout) :: ikpt10,npw1,npw2
2731  real(dp),intent(in) :: ecut1,ecut2,ecut2_eff
2732  type(MPI_type),intent(inout) :: mpi_enreg1,mpi_enreg2
2733 !arrays
2734  integer,intent(in) :: indkk(nkpt2*sppoldbl,6),istwfk1(nkpt1),istwfk2(nkpt2)
2735  integer,intent(in) :: ngfft1(18),ngfft2(18),symrel(3,3,nsym)
2736  integer,intent(inout) :: kg1(3,mpw1),kg2(3,mpw2)
2737  real(dp),intent(in) :: gmet1(3,3),gmet2(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2)
2738  real(dp),intent(in) :: rprimd2(3,3),tnons(3,nsym)
2739  real(dp),intent(inout) :: cg1(2,mcg1),cg2(2,mcg2)
2740  real(dp),intent(inout) :: eig_k1(mband1*(2*mband1)**formeig)
2741  real(dp),intent(inout) :: eig_k2(mband2*(2*mband2)**formeig),occ_k1(mband1)
2742  real(dp),intent(inout) :: occ_k2(mband2)
2743 
2744 !Local variables ------------------------------
2745 !scalars
2746  integer,parameter :: nkpt_max=50,tobox=1,tosph=-1
2747  integer :: conv_tnons,convert,fftalg,fold1,fold2,foldim,foldre,i1,i2,iband
2748  integer :: iband_first,iband_last,icgmod,ierr,index,ipw
2749  integer :: ispinor,ispinor1,ispinor2,ispinor_first,ispinor_last
2750  integer :: istwf10_k,istwf1_k,istwf2_k,isym,itimrev,jsign
2751  integer :: mgfft1,mgfft2,n1,n2,n3,n4,n5,n6
2752  integer :: nbremn,npwtot,nspinor_index,nspinor1_this_proc,nspinor2_this_proc
2753  integer :: order,ortalgo,seed
2754  real(dp) :: ai,ar,arg,bi,br,eig_tmp,spinrots,spinrotx,spinroty,spinrotz
2755  character(len=500) :: message
2756  integer, parameter :: int64 = selected_int_kind(18)
2757  !arrays
2758  integer :: atindx(1),identity(3,3),ngfft_now(18),no_shift(3),shiftg(3)
2759  integer :: symm(3,3),symrel_conv(3,3)
2760  integer,allocatable :: gbound1(:,:),gbound2(:,:)
2761  real(dp) :: kpoint1(3),kpoint2_sph(3),phktnons(2,1),spinrot(4),tnons_conv(3),tsec(2)
2762  real(dp),allocatable :: cfft(:,:,:,:),dum(:,:),phase1d(:,:),phase3d(:,:)
2763  real(dp),allocatable :: wavef1(:,:),wavef2(:,:),wavefspinor(:,:)
2764 
2765 ! *************************************************************************
2766 
2767  mgfft1=maxval(ngfft1(1:3))
2768  mgfft2=maxval(ngfft2(1:3))
2769  if(.false.)write(std_out,*)occ_k1 ! just to keep occ_k1 as an argument before resolving the issue of its transfer
2770 
2771  if(nspinor1/=1 .and. nspinor1/=2)then
2772    write(message,'(a,i0)')'The argument nspinor1 must be 1 or 2, while it is nspinor1 = ',nspinor1
2773    MSG_BUG(message)
2774  end if
2775 
2776  if(nspinor2/=1 .and. nspinor2/=2)then
2777    write(message,'(a,i0)')' The argument nspinor2 must be 1 or 2, while it is nspinor2=',nspinor2
2778    MSG_BUG(message)
2779  end if
2780 
2781  if(nspinor1==2 .and. mod(nbd1,2)/=0)then
2782    write(message,'(a,i0)')' When nspinor1 is 2, nbd1 must be even, while it is nbd1 = ',nbd1
2783    MSG_BUG(message)
2784  end if
2785 
2786  if(nspinor2==2 .and. mod(nbd2,2)/=0)then
2787    write(message,'(a,i0)')'  When nspinor2 is 2, nbd2 must be even, while it is nbd2=',nbd2
2788    MSG_BUG(message)
2789  end if
2790 
2791  if(nbd1/nspinor1>nbd2/nspinor2)then
2792    write(message, '(3a,2i6,3a,2i6,a)' )&
2793 &   'In wfconv, the nbd/nspinor ratio cannot decrease. However,',ch10,&
2794 &   'the initial quantities are nbd1,nspinor1=',nbd1,nspinor1,', and',ch10,&
2795 &   'the requested final quantities are nbd2,nspinor2=',nbd2,nspinor2,'.'
2796    MSG_BUG(message)
2797  end if
2798 
2799  ngfft_now(1:3)=ngfft1(1:3)
2800  ngfft_now(8:18)=ngfft1(8:18)
2801 !This line is the reason why ngfft_now has to be introduced
2802  ngfft_now(7)=101
2803  ngfft_now(4:6)=ngfft_now(1:3)
2804  n1=ngfft_now(1) ; n2=ngfft_now(2) ; n3=ngfft_now(3)
2805  n4=ngfft_now(4) ; n5=ngfft_now(5) ; n6=ngfft_now(6)
2806  fftalg=ngfft_now(7)
2807 
2808 !Parallelization over spinors management
2809  nspinor1_this_proc=max(1,nspinor1/mpi_enreg1%nproc_spinor)
2810  nspinor2_this_proc=max(1,nspinor2/mpi_enreg2%nproc_spinor)
2811 
2812 !In order to generate IN PLACE new wfs from old wfs, the loop
2813 !over bands and spinors must be done in one direction or the other,
2814 !depending on npw1 and npw2, nspinor1 and nspinor2.
2815 !If nspinor1=1 and nspinor2=2 , note that one will generate
2816 !from nbd1 states of npw1 coefficients,
2817 !2*nbd1 states of 2*npw2 coefficients. nbd1 cancels in comparing
2818 !these expressions, but nspinor2 appears squared.
2819 !The same line of thought works for the case nspinor1=2 and nspinor2=1
2820  order=1
2821  iband_first=1   ; iband_last=nbd1
2822  ispinor_first=1 ; ispinor_last=nspinor1
2823  if(nspinor1==2 .and. nspinor2==1)then
2824    order=2 ; iband_last=nbd1-1 ; ispinor_last=1
2825  end if
2826 !Here, reverse the order if needed
2827  if( npw2*nspinor2**2 > npw1*nspinor1**2 )then
2828    order=-order
2829    iband_first=iband_last       ; iband_last=1
2830    ispinor_first=ispinor_last   ; ispinor_last=1
2831  end if
2832 
2833  kpoint1(:)=kptns1(:,ikpt1)
2834  istwf1_k=istwfk1(ikpt1)
2835 
2836  kpoint2_sph(:)=0.0_dp
2837  if(ceksp2==0)kpoint2_sph(:)=kptns2(:,ikpt2)
2838  istwf2_k=istwfk2(ikpt2)
2839 
2840 !DEBUG
2841 !write(std_out,*)'ecut1,ecut2_eff=',ecut1,ecut2_eff
2842 !write(std_out,*)'gmet1,gmet2=',gmet1,gmet2
2843 !write(std_out,*)'kpoint1,kpoint2_sph=',kpoint1,kpoint2_sph
2844 !write(std_out,*)'nspinor1,nspinor2',nspinor1,nspinor2
2845 !write(std_out,*)'istwf1_k,istwf2_k=',istwf1_k,istwf2_k
2846 !write(std_out,*)'nbd1,tol8=',nbd1,tol8
2847 !ENDDEBUG
2848 
2849 !Determine whether it will be needed to convert the existing
2850 !wavefunctions, or simply to complete them.
2851 
2852  convert=0
2853  if(nbd1/=0)then
2854    if(abs(ecut2_eff-ecut1)>tol8)convert=convert+1
2855    if(sum(abs(gmet2(:,:)-gmet1(:,:)))>tol8)convert=convert+2
2856    if(sum(abs(kpoint2_sph(:)-kpoint1(:)))>tol8)convert=convert+4
2857    if(nspinor2/=nspinor1)convert=convert+8
2858    if(istwf2_k/=istwf1_k)convert=convert+16
2859  end if
2860 
2861 !This is a supplementary check
2862  if(restart==1 .and. convert/=0)then
2863    MSG_BUG('Restart==1 and convert/=0 are exclusive')
2864  end if
2865 
2866 !Determine whether symmetries must be used
2867  conv_tnons=0
2868  no_shift(:)=0
2869  identity(:,:)=0
2870  identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
2871  isym=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,2)
2872 !DEBUG
2873 !write(std_out,*)' wfconv : isym=',isym
2874 !ENDDEBUG
2875  itimrev=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,6)
2876  if(isym/=0)then
2877    symrel_conv(:,:)=symrel(:,:,isym)
2878    call mati3inv(symrel_conv,symm)
2879    shiftg(:)=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,3:5)
2880    tnons_conv(:)=tnons(:,isym)
2881    if(sum(tnons_conv(:)**2)>tol8)then
2882 !    Need to compute phase factors associated with nonsymmorphic translations.
2883      conv_tnons=1
2884      ABI_ALLOCATE(phase3d,(2,npw1))
2885      ABI_ALLOCATE(phase1d,(2,(2*n1+1)+(2*n2+1)+(2*n3+1)))
2886 !    Although the routine getph is originally written for
2887 !    atomic phase factors, it does precisely what we want
2888      atindx(1)=1
2889      call getph(atindx,1,n1,n2,n3,phase1d,tnons_conv)
2890    end if
2891    if(nspinor1==2 .and. nspinor2==2)then
2892 !    Compute rotation in spinor space
2893      call getspinrot(rprimd2,spinrot,symrel_conv)
2894    end if
2895  else
2896    shiftg(:)=0
2897    symm(:,:)=identity(:,:)
2898    spinrot(:)=zero
2899    spinrot(1)=one
2900  end if
2901  if(itimrev/=0)then
2902    symm(:,:)=-symm(:,:)
2903  end if
2904 
2905 !DEBUG
2906 !write(std_out,'(a,i3,2x,3i3,2x,9i3)')' wfconv : isym,shiftg,symm=',isym,shiftg,symm
2907 !write(std_out,*)' wfconv : ecut2_eff,ecut1=',ecut2_eff,ecut1
2908 !write(std_out,*)' wfconv : istwf1_k,istwf2_k=',istwf1_k,istwf2_k
2909 !write(std_out,*)' wfconv : kpoint1(:),kpoint2_sph(:)=',&
2910 !& kpoint1(:),kpoint2_sph(:)
2911 !write(std_out,*)' wfconv : nspinor1,nspinor2=',nspinor1,nspinor2
2912 !ENDDEBUG
2913 
2914 !if (mpi_enreg1%fft_option_lob==0) mpi_enreg1%fft_option_lob=1
2915 !if (mpi_enreg2%fft_option_lob==0) mpi_enreg2%fft_option_lob=1
2916 
2917  if (restart==2.and.(convert/=0.or.(nbd2/nspinor2>nbd1/nspinor1.and.formeig==0))) then
2918 !  kg2 is needed both for FFT grid conversion and for envlop
2919 !  Choose the center of the sphere : either gamma, or each k-point
2920    kpoint2_sph(:)=0.0_dp
2921    if(ceksp2==0)kpoint2_sph(:)=kptns2(:,ikpt2)
2922    istwf2_k=istwfk2(ikpt2)
2923    call kpgsph(ecut2_eff,exchn2n3d,gmet2,0,ikpt2,istwf2_k,kg2,kpoint2_sph,1,mpi_enreg2,mpw2,npw2)
2924  end if
2925 
2926  if(convert/=0)then
2927    istwf10_k=0
2928    if(ikpt10/=0)istwf10_k=istwfk1(ikpt10)
2929 
2930 !  Only need G sphere if different from last time
2931    if ( ikpt1/=ikpt10 .or. istwf1_k/=istwf10_k ) then
2932 
2933      call kpgsph (ecut1,exchn2n3d,gmet1,0,ikpt1,istwf1_k,kg1,kpoint1,1,mpi_enreg1,mpw1,npw1)
2934      if (debug>0) then
2935        write(message, '(a,f8.3,a,a,3f8.5,a,a,i3,a,3(a,3es16.8,a),a,3i4,a,i5,a)' )&
2936 &       ' wfconv: called kpgsph with ecut1=',ecut1,ch10,&
2937 &       '  kpt1=',kptns1(1:3,ikpt1),ch10,&
2938 &       '  istwf1_k=',istwf1_k,ch10,&
2939 &       '  gmet1= ',gmet1(1:3,1),ch10,&
2940 &       '         ',gmet1(1:3,2),ch10,&
2941 &       '         ',gmet1(1:3,3),ch10,&
2942 &       '  ngfft=',ngfft_now(1:3),' giving npw1=',npw1,'.'
2943        call wrtout(std_out,message,'PERS')
2944      end if
2945      ikpt10 = ikpt1
2946      istwf10_k=istwf1_k
2947    end if
2948 
2949    if(conv_tnons==1)then
2950      arg=two_pi*(kpoint1(1)*tnons_conv(1)+ kpoint1(2)*tnons_conv(2)+ kpoint1(3)*tnons_conv(3) )
2951      phktnons(1,1)=cos(arg)
2952      phktnons(2,1)=sin(arg)
2953 !    Convert 1D phase factors to 3D phase factors exp(i 2 pi (k+G).tnons )
2954      call ph1d3d(1,1,kg1,1,1,npw1,n1,n2,n3,phktnons,phase1d,phase3d)
2955    end if
2956 
2957    ABI_ALLOCATE(cfft,(2,n4,n5,n6))
2958    ABI_ALLOCATE(wavef1,(2,npw1))
2959    ABI_ALLOCATE(wavef2,(2,npw2))
2960    if(nspinor1==2 .and. nspinor2==2) then
2961      ABI_ALLOCATE(wavefspinor,(2,2*npw2))
2962    end if
2963    ABI_ALLOCATE(gbound1,(2*mgfft1+8,2))
2964    ABI_ALLOCATE(gbound2,(2*mgfft2+8,2))
2965    call sphereboundary(gbound1,istwf1_k,kg1,mgfft1,npw1)
2966    call sphereboundary(gbound2,istwf2_k,kg2,mgfft2,npw2)
2967 
2968 !  Take old wf from sphere->box, the new from box->sphere
2969 !  One pays attention not to have a problem of erasing data when replacing
2970 !  a small set of coefficient by a large set, or the reverse.
2971 !  This is the reason of the use of order, _first and _last variables,
2972 !  defined earlier.
2973    nspinor_index=mpi_enreg1%me_spinor+1
2974    do iband=iband_first,iband_last,order
2975      do ispinor1=ispinor_first,ispinor_last,order
2976        ispinor=ispinor1
2977        if (mpi_enreg1%paral_spinor==1) then
2978          if (ispinor1==nspinor_index) then
2979            ispinor=1
2980          else
2981            if (nspinor1==2.and.nspinor2==2) wavefspinor(:,(ispinor1-1)*npw2+1:ispinor1*npw2)=zero
2982            cycle
2983          end if
2984        end if
2985 
2986 !      Copy input wf
2987        i1=(ispinor-1)*npw1+(iband-1)*nspinor1_this_proc*npw1+icg1
2988        wavef1(:,1:npw1)=cg1(:,i1+1:i1+npw1)
2989 
2990 !      Make symmetry-induced conversion, if needed (translation part)
2991        if(conv_tnons==1)then
2992 !$OMP PARALLEL DO PRIVATE(ai,ar)
2993          do ipw=1,npw1
2994            ar=phase3d(1,ipw)*wavef1(1,ipw)-phase3d(2,ipw)*wavef1(2,ipw)
2995            ai=phase3d(2,ipw)*wavef1(1,ipw)+phase3d(1,ipw)*wavef1(2,ipw)
2996            wavef1(1,ipw)=ar
2997            wavef1(2,ipw)=ai
2998          end do
2999        end if
3000 
3001 !      Take into account time-reversal symmetry, if needed, in the scalar case
3002        if(itimrev==1 .and. (nspinor1==1 .or. nspinor2==1))then
3003 !$OMP PARALLEL DO
3004          do ipw=1,npw1
3005            wavef1(2,ipw)=-wavef1(2,ipw)
3006          end do
3007        end if
3008 
3009 !      DEBUG
3010 !      write(std_out,*)' wfconv : before sphere, isym,ispinor=',isym,ispinor
3011 !      write(std_out,*)' no_shift,identity=',no_shift,identity
3012 !      write(std_out,*)' shiftg,symm=',shiftg,symm
3013 !      stop
3014 !      This debugging sequence is an attempt to rotate spinors,
3015 !      and works indeed for test13, when symmetry 9 is used ...
3016 !      if(isym==9 .and. ispinor==1)then
3017 !      write(std_out,*)' wfconv : gives a 120 degree rotation to first component'
3018 !      do ipw=1,npw1
3019 !      ar=-            half*wavef1(1,ipw)-sqrt(three)*half*wavef1(2,ipw)
3020 !      ai= sqrt(three)*half*wavef1(1,ipw)-            half*wavef1(2,ipw)
3021 !      wavef1(1,ipw)=ar
3022 !      wavef1(2,ipw)=ai
3023 !      end do
3024 !      end if
3025 !      ENDDEBUG
3026 
3027 !      Convert wf, and also include the symmetry operation and shiftg.
3028        call sphere(wavef1,1,npw1,cfft,n1,n2,n3,n4,n5,n6,kg1,istwf1_k,tobox,&
3029 &       mpi_enreg1%me_g0,no_shift,identity,one)
3030 
3031        call sphere(wavef2,1,npw2,cfft,n1,n2,n3,n4,n5,n6,kg2,istwf2_k,tosph,&
3032 &       mpi_enreg2%me_g0,shiftg,symm,one)
3033 
3034        if(nspinor2==1 )then
3035          i2=(ispinor-1)*npw2+(iband-1)*nspinor2_this_proc*npw2+icg2
3036          cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2)
3037        else if(nspinor1==2.and.nspinor2==2)then
3038 !        Will treat this case outside of the ispinor loop
3039          i2=(ispinor1-1)*npw2
3040          wavefspinor(:,i2+1:i2+npw2)=wavef2(:,1:npw2)
3041        else if(nspinor1==1 .and. nspinor2==2)then
3042 !        The number of bands is doubled, and the number of coefficients
3043 !        is doubled also
3044          if (mpi_enreg2%paral_spinor==0) then
3045            i2=(iband-1)*nspinor2_this_proc*nspinor2_this_proc*npw2+icg2
3046            cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2)
3047            cg2(:,i2+npw2+1:i2+2*npw2)=zero
3048            cg2(:,i2+2*npw2+1:i2+3*npw2)=zero
3049            cg2(:,i2+3*npw2+1:i2+4*npw2)=wavef2(:,1:npw2)
3050          else
3051            i2=(iband-1)*nspinor2_this_proc*npw2+icg2
3052            if (nspinor_index==1) then
3053              cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2)
3054              cg2(:,i2+npw2+1:i2+2*npw2)=zero
3055            else
3056              cg2(:,i2+1:i2+npw2)=zero
3057              cg2(:,i2+npw2+1:i2+2*npw2)=wavef2(:,1:npw2)
3058            end if
3059          end if
3060        end if
3061      end do ! ispinor=ispinor_first,ispinor_last,order
3062 
3063      if(nspinor1==2.and.nspinor2==2)then
3064 !      Take care of possible parallelization over spinors
3065        if (mpi_enreg2%paral_spinor==1) then
3066          call xmpi_sum(wavefspinor,mpi_enreg2%comm_spinor,ierr)
3067        end if
3068 !      Take care of time-reversal symmetry, if needed
3069        if(itimrev==1)then
3070 !        Exchange spin-up and spin-down
3071 !        Make complex conjugate of one component,
3072 !        and change sign of other component
3073 !$OMP PARALLEL DO PRIVATE(ipw,ar,ai) SHARED(wavefspinor,npw2)
3074          do ipw=1,npw2
3075 !          Here, change sign of real part
3076            ar=-wavefspinor(1,ipw)
3077            ai= wavefspinor(2,ipw)
3078            wavefspinor(1,ipw)= wavefspinor(1,npw2+ipw)
3079 !          Here, change sign of imaginary part
3080            wavefspinor(2,ipw)=-wavefspinor(2,npw2+ipw)
3081            wavefspinor(1,npw2+ipw)=ar
3082            wavefspinor(2,npw2+ipw)=ai
3083          end do
3084        end if ! itimrev==1
3085 
3086 !      Rotation in spinor space
3087 !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(npw2,spinrot,wavefspinor)
3088        spinrots=spinrot(1)
3089        spinrotx=spinrot(2)
3090        spinroty=spinrot(3)
3091        spinrotz=spinrot(4)
3092 !$OMP DO
3093        do ipw=1,npw2
3094          ar=wavefspinor(1,ipw)
3095          ai=wavefspinor(2,ipw)
3096          br=wavefspinor(1,npw2+ipw)
3097          bi=wavefspinor(2,npw2+ipw)
3098          wavefspinor(1,ipw)     = spinrots*ar-spinrotz*ai +spinroty*br-spinrotx*bi
3099          wavefspinor(2,ipw)     = spinrots*ai+spinrotz*ar +spinroty*bi+spinrotx*br
3100          wavefspinor(1,npw2+ipw)=-spinroty*ar-spinrotx*ai +spinrots*br+spinrotz*bi
3101          wavefspinor(2,npw2+ipw)=-spinroty*ai+spinrotx*ar +spinrots*bi-spinrotz*br
3102        end do
3103 !$OMP END DO
3104 !$OMP END PARALLEL
3105 
3106 !      Save wavefunction
3107        i2=(iband-1)*nspinor2_this_proc*npw2+icg2
3108        if (mpi_enreg2%paral_spinor==0) then
3109          cg2(:,i2     +1:i2+  npw2)=wavefspinor(:,1:npw2)
3110          cg2(:,i2+npw2+1:i2+2*npw2)=wavefspinor(:,npw2+1:2*npw2)
3111        else
3112          if (nspinor_index==1) then
3113            cg2(:,i2+1:i2+npw2)=wavefspinor(:,1:npw2)
3114          else
3115            cg2(:,i2+1:i2+npw2)=wavefspinor(:,npw2+1:2*npw2)
3116          end if
3117        end if
3118      end if ! nspinor1==2 .and. nspinor2==2
3119 
3120    end do
3121 
3122 !  Take care of copying eig and occ when nspinor increases or decreases
3123    if(nspinor1==1.and.nspinor2==2)then
3124      if(formeig==0)then
3125 !      Note the reverse order, needed in case inplace=1
3126        do iband=nbd1,1,-1
3127 !        use eig_tmp to avoid bug on ifort10.1 x86_64
3128          eig_tmp=eig_k1(iband)
3129          eig_k2(2*iband-1:2*iband)=eig_tmp
3130 !        occ_tmp=occ_k1(iband)*0.5_dp
3131 !        occ_k2(2*iband-1:2*iband )=occ_tmp
3132        end do
3133      else
3134        call wrtout(std_out,' wfconv: not yet coded, formeig=1!',"COLL")
3135      end if
3136    end if
3137    if(nspinor1==2 .and. nspinor2==1)then
3138      if(formeig==0)then
3139        do iband=1,nbd1
3140 !        use eig_tmp to avoid bug on ifort10.1 x86_64
3141          eig_tmp=eig_k1(2*iband-1)
3142          eig_k2(iband)=eig_tmp
3143 !        occ_tmp=occ_k1(2*iband-1)*2.0_dp
3144 !        occ_k2(iband)=occ_tmp
3145        end do
3146      else
3147        call wrtout(std_out,' wfconv: not yet coded, formeig=1!',"COLL")
3148      end if
3149    end if
3150 
3151    ABI_DEALLOCATE(cfft)
3152    ABI_DEALLOCATE(gbound1)
3153    ABI_DEALLOCATE(gbound2)
3154    ABI_DEALLOCATE(wavef1)
3155    ABI_DEALLOCATE(wavef2)
3156    if(nspinor1==2 .and. nspinor2==2) then
3157      ABI_DEALLOCATE(wavefspinor)
3158    end if
3159 
3160  else if(convert==0)then
3161 
3162    if(inplace==0)then
3163 !    Must copy cg, eig and occ if not in-place while convert==0
3164 !    Note that npw1=npw2, nspinor1=nspinor2
3165      cg2(:,1+icg2:npw1*nspinor1_this_proc*nbd1+icg2)=&
3166 &     cg1(:,1+icg1:npw1*nspinor1_this_proc*nbd1+icg1)
3167      eig_k2(:)=eig_k1(:)
3168 !    occ_k2(:)=occ_k1(:)
3169    end if
3170 
3171  end if ! End of if convert/=0
3172 
3173  if(conv_tnons==1) then
3174    ABI_DEALLOCATE(phase1d)
3175    ABI_DEALLOCATE(phase3d)
3176  end if
3177 
3178 
3179 !If not enough bands, complete with random numbers or zeros
3180  if(nbd2/nspinor2>nbd1/nspinor1)then
3181    if(formeig==0)then
3182 
3183 !    Ground state wf and eig case
3184      eig_k2((nbd1/nspinor1)*nspinor2+1:nbd2)=huge(0.0_dp)/10.0_dp
3185      occ_k2((nbd1/nspinor1)*nspinor2+1:nbd2)=0.0_dp
3186      index=(nbd1/nspinor1)*nspinor2*npw2*nspinor2_this_proc
3187 
3188 !    Initialisation of wavefunctions
3189 !    One needs to initialize wfs in such a way to avoid symmetry traps,
3190 !    and to avoid linear dependencies between wavefunctions
3191 !    No need for a difference for different k points and/or spin-polarization
3192 
3193      if (mpi_enreg1%paral_kgb == 1) then
3194        npwtot=npw2
3195        call timab(539,1,tsec)
3196        call xmpi_sum(npwtot, mpi_enreg1%comm_bandfft, ierr)
3197        call timab(539,2,tsec)
3198      end if
3199 
3200      do iband=(nbd1/nspinor1)*nspinor2+1,nbd2
3201        do ispinor2=1,nspinor2_this_proc
3202          ispinor=ispinor2;if (nspinor2_this_proc/=nspinor2) ispinor=mpi_enreg2%me_spinor+1
3203          jsign=1;if (ispinor==2) jsign=-1
3204 
3205          do ipw=1,npw2
3206            index=index+1
3207 !          Different seed for different planewave and band
3208 !          DEBUG seq==par
3209 !          if(.false.) then
3210 !          ENDDEBUG seq==par
3211 
3212            if ( mpi_enreg2%paral_kgb /= 1) then
3213              seed=(iband-1)*npw2*nspinor2 + (ispinor-1)*npw2 + ipw
3214            else
3215              seed=jsign*(iband*(kg2(1,ipw)*npwtot*npwtot + kg2(2,ipw)*npwtot + kg2(3,ipw)))
3216            end if
3217 
3218            if(randalg == 0) then
3219 !            For portability, use only integer numbers
3220 !            The series of couples (fold1,fold2) is periodic with a period of
3221 !            3x5x7x11x13x17x19x23x29x31, that is, larger than 2**32, the largest integer*4
3222 !            fold1 is between 0 and 34, fold2 is between 0 and 114. As sums of five
3223 !            uniform random variables, their distribution is close to a gaussian
3224              fold1=mod(seed,3)+mod(seed,5)+mod(seed,7)+mod(seed,11)+mod(seed,13)
3225              fold2=mod(seed,17)+mod(seed,19)+mod(seed,23)+mod(seed,29)+mod(seed,31)
3226 !            The gaussian distributions are folded, in order to be back to a uniform distribution
3227 !            foldre is between 0 and 20, foldim is between 0 and 18
3228              foldre=mod(fold1+fold2,21)
3229              foldim=mod(3*fold1+2*fold2,19)
3230              cg2(1,index+icg2)=dble(foldre)
3231              cg2(2,index+icg2)=dble(foldim)
3232            else
3233              ! (AL) Simple linear congruential generator from
3234              ! numerical recipes, modulo'ed and 64bit'ed to avoid
3235              ! overflows (NAG doesn't like overflows, even though
3236              ! they are perfectly legitimate here). Then, we get some
3237              ! lowest order bits and sum them, as the previous
3238              ! generator, to get quasi-normal numbers.
3239              ! This is clearly suboptimal and might cause problems,
3240              ! but at least it doesn't seem to create linear
3241              ! dependencies and local minima like the previous one.
3242              ! it's not trivial to generate good reproductible random
3243              ! numbers in parallel. Patches welcome !
3244              ! Note a fun fortran fact : MOD simply ignores 64 bits integer
3245              ! and casts them into 32bits, so we use MODULO.
3246              fold1 = modulo(1664525_int64 * seed  + 1013904223_int64, 2147483648_int64)
3247              fold2 = modulo(1664525_int64 * fold1 + 1013904223_int64, 2147483648_int64)
3248              fold1=modulo(fold1,3)+modulo(fold1,5)+modulo(fold1,7)+modulo(fold1,11)+modulo(fold1,13)
3249              fold2=modulo(fold2,3)+modulo(fold2,5)+modulo(fold2,7)+modulo(fold2,11)+modulo(fold2,13)
3250              cg2(1,index+icg2)=dble(fold1)/34-0.5
3251              cg2(2,index+icg2)=dble(fold2)/34-0.5
3252            end if
3253          end do
3254        end do
3255 
3256 !      XG030513 : MPIWF need to impose cg to zero when at Gamma
3257 !      Time-reversal symmetry for k=gamma impose zero imaginary part at G=0
3258 !      XG : I do not know what happens for spin-orbit here :
3259        if(istwf2_k==2 .and. mpi_enreg2%me_g0==1) then
3260          cg2(2,1+(iband-1)*npw2*nspinor2_this_proc+icg2)=zero
3261        end if
3262      end do
3263 
3264 !    Multiply with envelope function to reduce kinetic energy
3265      icgmod=icg2+npw2*nspinor2_this_proc*(nbd1/nspinor1)
3266      nbremn=nbd2-nbd1
3267      call cg_envlop(cg2,ecut2,gmet2,icgmod,kg2,kpoint2_sph,mcg2,nbremn,npw2,nspinor2_this_proc)
3268 
3269      if(ikpt2<=nkpt_max)then
3270        write(message,'(3(a,i6))')' wfconv:',nbremn,' bands initialized randomly with npw=',npw2,', for ikpt=',ikpt2
3271        call wrtout(std_out,message,'PERS')
3272      end if
3273 
3274    else if(formeig==1)then
3275 
3276 !    For response function, put large numbers in the remaining of the
3277 !    eigenvalue array (part of it was already filled in calling routine)
3278 !    WARNING : Change of nspinor not yet coded
3279      eig_k2(1+2*nbd1*nbd2 : 2*nbd2*nbd2)=huge(0.0_dp)/10.0_dp
3280 !    Initialisation of wfs with 0 s
3281      index=npw2*nbd1*nspinor2_this_proc
3282      do iband=nbd1+1,nbd2
3283        do ipw=1,npw2*nspinor2_this_proc
3284          index=index+1
3285          cg2(:,index+icg2)=zero
3286        end do
3287      end do
3288 
3289      if(ikpt2<=nkpt_max)then
3290        nbremn=nbd2-nbd1
3291        write(message,'(a,i6,a,i7,a,i4)')' wfconv :',nbremn,' bands set=0 with npw=',npw2,', for ikpt=',ikpt2
3292        call wrtout(std_out,message,'PERS')
3293      end if
3294 
3295    end if ! End of initialisation to 0
3296  end if
3297 
3298 !Orthogonalize GS wfs
3299  !if (.False.) then
3300  if (optorth==1.and.formeig==0.and.mpi_enreg2%paral_kgb/=1) then
3301    ABI_ALLOCATE(dum,(2,0))
3302    ortalgo=0 !;ortalgo=3
3303    call pw_orthon(icg2,0,istwf2_k,mcg2,0,npw2*nspinor2_this_proc,nbd2,ortalgo,dum,0,cg2,&
3304 &   mpi_enreg2%me_g0,mpi_enreg2%comm_bandspinorfft)
3305    ABI_DEALLOCATE(dum)
3306  end if
3307 
3308 end subroutine wfconv

m_inwffil/initwf [ Functions ]

[ Top ] [ m_inwffil ] [ Functions ]

NAME

 initwf

FUNCTION

 Initialization of wavefunctions.
 If formeig==1, and partially filled case, I am not sure that the eig_k are initialized properly ...
 formeig option (format of the eigenvalues and eigenvector) :
   0 => ground-state format (initialisation of
        eigenvectors with random numbers, vector of eigenvalues)
   1 => respfn format (initialisation of
        eigenvectors with 0 s, hermitian matrix of eigenvalues)

INPUTS

 formeig=see above
 headform=header format (might be needed to read the block of wfs)
 icg=shift to be given to the location of the data in the array cg
 ikpt= number of the k point of which the wf is initialised
 spin=spin index
 mcg=dimension of the cg array
 mpi_enreg=information about MPI parallelization
 nband_k=number of bands at this particular k point
 nkpt=number of k points
 npw=number of plane waves
 nspinor=number of spinorial components of the wavefunctions (on current proc)
 wff1=structure info for file containing wavefunctions (when needed)

OUTPUT

 cg(2,mcg)=complex wf array
  if ground state format (formeig=0):
    eig_k(nband_k)=list of eigenvalues (input or init to large number), hartree
  if respfn format (formeig=1):
    eig_k(2*nband_k*nband_k)= matrix of eigenvalues (input or init to large number), hartree

SIDE EFFECTS

 Input/output:
 occ_k(nband_k)=list of occupations (input or left to their initial value)
 ikptsp_old=number of the previous spin-k point, or 0 if first call of present file

PARENTS

      wfsinp

CHILDREN

      rwwf,timab,wffreadskipk,wfk_close,wfk_open_read,wfk_read_band_block
      wrtout

SOURCE

1863 subroutine initwf(cg,eig_k,formeig,headform,icg,ikpt,ikptsp_old,&
1864 &  spin,mcg,mpi_enreg,nband_k,nkpt,npw,nspinor,occ_k,wff1)
1865 
1866 
1867 !This section has been created automatically by the script Abilint (TD).
1868 !Do not modify the following lines by hand.
1869 #undef ABI_FUNC
1870 #define ABI_FUNC 'initwf'
1871 !End of the abilint section
1872 
1873  implicit none
1874 
1875 !Arguments ------------------------------------
1876 !scalars
1877  integer,intent(in) :: formeig,headform,icg,ikpt,spin,mcg,nband_k,nkpt,npw,nspinor
1878  integer,intent(inout) :: ikptsp_old
1879  type(MPI_type),intent(in) :: mpi_enreg
1880  type(wffile_type),intent(inout) :: wff1
1881 !arrays
1882  real(dp),intent(inout) :: occ_k(nband_k)
1883  real(dp),intent(inout) :: cg(2,mcg),eig_k((2*nband_k)**formeig*nband_k) !vz_i
1884 
1885 !Local variables-------------------------------
1886 !scalars
1887  integer,parameter :: nkpt_max=50
1888  integer :: ikpt0,nband_disk,tim_rwwf
1889  character(len=500) :: msg
1890 !arrays
1891  integer,allocatable :: kg_dum(:,:)
1892  real(dp) :: tsec(2)
1893 #if 0
1894  integer :: iomode,comm,funt,ierr
1895  type(wfk_t) :: Wfk
1896 #endif
1897 
1898 ! *************************************************************************
1899 
1900 !DEBUG
1901 !write(std_out,*)' initwf : enter, ikptsp_old,ikpt,spin,nkpt= ',ikptsp_old,ikpt,spin,nkpt
1902 !stop
1903 !ENDDEBUG
1904 
1905 #if 0
1906  MSG_WARNING("Entering new IO section")
1907 !call WffClose(wff1,ierr)
1908  comm   = MPI_enreg%comm_cell
1909  iomode = iomode_from_fname(wff1%fname)
1910  call wfk_open_read(Wfk,wff1%fname,formeig,iomode,get_unit(),comm)
1911  call wfk_read_band_block(Wfk,(/1,nband_k/),ikpt,spin,xmpio_at,cg_k=cg(1:,icg:),eig_k=eig_k,occ_k=occ_k)
1912  call wfk_close(Wfk)
1913 !call clsopn(wff1)
1914  RETURN
1915 #endif
1916 
1917  call timab(770,1,tsec)
1918  call timab(771,1,tsec)
1919 
1920  ABI_ALLOCATE(kg_dum,(3,0))
1921 
1922 !Skip wavefunctions for k-points not treated by this proc.
1923 !(from ikptsp_old+1 to ikpt+(spin-1)*nkpt-1)
1924  if (ikptsp_old<ikpt+(spin-1)*nkpt-1) then
1925    do ikpt0=ikptsp_old+1,ikpt+(spin-1)*nkpt-1
1926      call WffReadSkipK(formeig,headform,ikpt0,spin,mpi_enreg,wff1)
1927    end do
1928  end if
1929 
1930 !DEBUG
1931 !write(std_out,*)' initwf : before rwwf'
1932 !write(std_out,*)' formeig,icg,ikpt,spin=',formeig,icg,ikpt,spin
1933 !write(std_out,*)' nband_k,nband_disk,npw,nspinor=',nband_k,nband_disk,npw,nspinor
1934 !write(std_out,*)' unwff1=',unwff1
1935 !stop
1936 !ENDDEBUG
1937 
1938  if(mpi_enreg%paralbd==0)tim_rwwf=2
1939  if(mpi_enreg%paralbd==1)tim_rwwf=20
1940 
1941  call timab(771,2,tsec)
1942 
1943  call rwwf(cg,eig_k,formeig,headform,icg,ikpt,spin,kg_dum,nband_k,mcg,mpi_enreg,nband_k,nband_disk,&
1944 & npw,nspinor,occ_k,1,0,tim_rwwf,wff1)
1945 
1946  call timab(772,1,tsec)
1947 
1948  if (ikpt<=nkpt_max) then
1949    write(msg,'(3(a,i0))')' initwf: disk file gives npw= ',npw,' nband= ',nband_disk,' for kpt number= ',ikpt
1950    call wrtout(std_out,msg,'PERS')
1951  else if (ikpt==nkpt_max+1) then
1952    call wrtout(std_out,' initwf: the number of similar message is sufficient... stop printing them','PERS')
1953  end if
1954 
1955 !Check the number of bands on disk file against desired number. These are not required to agree)
1956  if (nband_disk/=nband_k) then
1957    write(msg,'(2(a,i0),3a,i0,3a)')&
1958 &   'For kpt number ',ikpt,' disk file has ',nband_disk,' bands',ch10,&
1959 &   'but input file gave nband= ',nband_k,'.',ch10,&
1960 &   'This is not fatal. Bands are skipped or filled with random numbers.'
1961    MSG_COMMENT(msg)
1962  end if
1963 
1964  if (ikpt<=nkpt_max) then
1965    write(msg,'(a,i0,a)')' initwf: ',nband_disk,' bands have been initialized from disk'
1966    call wrtout(std_out,msg,'PERS')
1967  end if
1968 
1969  ikptsp_old=ikpt+(spin-1)*nkpt
1970 
1971  ABI_DEALLOCATE(kg_dum)
1972 
1973  call timab(772,2,tsec)
1974  call timab(770,2,tsec)
1975 
1976 end subroutine initwf

m_inwffil/inwffil [ Functions ]

[ Top ] [ m_inwffil ] [ 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, ...)

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=information 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

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

m_inwffil/wfsinp [ Functions ]

[ Top ] [ m_inwffil ] [ Functions ]

NAME

 wfsinp

FUNCTION

 Do initialization of wavefunction files.
 Also call other relevant routines for this initialisation.
 Detailed description :
  - Initialize unit wff1 for input of wf data

 formeig option (format of the eigenvalues and occupations) :
   0 => ground-state format (initialisation of
        eigenvectors with random numbers, vector of eigenvalues,
        occupations are present)
   1 => respfn format (initialisation of
        eigenvectors with 0 s, hermitian matrix of eigenvalues)

INPUTS

  ecut0=kinetic energy cutoffs for basis sphere 0 (hartree) (if squeeze=1)
  ecut=kinetic energy cutoffs beyond which the coefficients of cg vanish (Ha)
   (needed only if squeeze=1)
  ecut_eff=effective kinetic energy planewave cutoff (hartree), needed
    to generate the sphere of plane wave
  exchn2n3d=if 1, n2 and n3 are exchanged
  formeig=explained above
  gmet(3,3), gmet0(3,3)=reciprocal space metrics (bohr^-2)
  headform0=header format (might be needed to read the block of wfs)
  indkk(nkpt*sppoldbl,6)=describe k point number of kptns0 that allows to
   generate wavefunctions closest to given kpt
   indkk(:,1)=k point number of kptns0
   indkk(:,2)=symmetry operation to be applied to kpt0, to give kpt0a
    (if 0, means no symmetry operation, equivalent to identity )
   indkk(:,3:5)=shift in reciprocal space to be given to kpt0a,
    to give kpt0b, that is the closest to kpt.
   indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
  indkk0(nkpt0,nkassoc)=list of k points that will be generated by k point number ikpt0
  istwfk(nkpt)=input parameter that describes the storage of wfs
  istwfk0(nkpt0)=input parameter that describes the storage of wfs in set0
  kptns(3,nkpt),kptns0(3,nkpt0)=k point sets (reduced coordinates)
  localrdwf=(for parallel case) if 1, the wff1%unwff file is local to each machine
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*nsppol
  mpi_enreg=information about MPI parallelization
  mpi_enreg0=information about MPI parallelization in set0
  mpw=maximum number of planewaves as dimensioned in calling routine
  mpw0=maximum number of planewaves on disk file
  nban_dp_rd(nkpt0*nsppol0)=number of bands to be read at each k point
  nband(nkpt*nsppol)=number of bands at each k point
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/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
  randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility
  restart= if 2, want conversion between wavefunctions
           if 1, direct restart is allowed (see hdr_check.f)
  rprimd(3,3)=dimensional primitive translations for real space (bohr)
  sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn
            if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries
  squeeze=1 if cg_disk is to be used, even with mkmem/=0
  symrel(3,3,nsym)=symmetry operations in real space in terms
   of primitive translations
  tnons(3,nsym)=nonsymmorphic translations for symmetry operations
  wff1, structure information for input and output files

OUTPUT

  if ground state format (formeig=0):
    eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
  if respfn format (formeig=1):
    eigen(2*mband*mband*nkpt*nsppol)=
         matrix of eigenvalues (input or init to large number), (Ha)
 Conditional output:
    cg_disk(2,mpw*nspinor*mband*mkmem*nsppol)=complex wf array
      be careful : an array of size cg(2,npw*nspinor), as used
      in the response function code, is not enough !

SIDE EFFECTS

  if ground state format (formeig=0):
    occ(mband*nkpt*nsppol)=occupations (from disk or left at their initial value)
    NOT OUTPUT NOW !

NOTES

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

WARNINGS

 For parallelism : no distinction yet between nban_dp_rd and nband

TODO

 THE DESCRIPTION IS TO BE COMPLETELY REVISED, AS THIS ONE COMES FROM inwffil.f

PARENTS

      inwffil

CHILDREN

      initwf,mpi_bcast,mpi_recv,mpi_send,pareigocc,timab,wfconv,wffreadskipk
      wrtout

SOURCE

1128 subroutine wfsinp(cg,cg_disk,ecut,ecut0,ecut_eff,eigen,exchn2n3d,&
1129 &                  formeig,gmet,gmet0,headform0,indkk,indkk0,istwfk,&
1130 &                  istwfk0,kptns,kptns0,localrdwf,mband,&
1131 &                  mcg,mcg_disk,mpi_enreg,mpi_enreg0,mpw,mpw0,nband,nban_dp_rd,&
1132 &                  ngfft,nkassoc,nkpt,nkpt0,npwarr,npwarr0,nspinor,&
1133 &                  nspinor0,nsppol,nsppol0,nsym,occ,optorth,prtvol,randalg,restart,rprimd,&
1134 &                  sppoldbl,squeeze,symrel,tnons,wff1)
1135 
1136 
1137 !This section has been created automatically by the script Abilint (TD).
1138 !Do not modify the following lines by hand.
1139 #undef ABI_FUNC
1140 #define ABI_FUNC 'wfsinp'
1141 !End of the abilint section
1142 
1143  implicit none
1144 
1145 !Arguments ------------------------------------
1146  integer, intent(in) :: exchn2n3d,formeig,headform0,localrdwf,mband,mcg,mcg_disk
1147  integer, intent(in) :: mpw,mpw0,nkassoc,nkpt,nkpt0,nspinor,nspinor0,nsppol,nsppol0,nsym
1148  integer, intent(in) :: optorth,prtvol,randalg,restart,sppoldbl,squeeze
1149  real(dp), intent(in) :: ecut,ecut0,ecut_eff
1150  type(MPI_type), intent(inout) :: mpi_enreg,mpi_enreg0
1151  type(wffile_type), intent(inout) :: wff1
1152  integer, intent(in) :: indkk(nkpt*sppoldbl,6),indkk0(nkpt0,nkassoc),istwfk(nkpt)
1153  integer, intent(in) :: istwfk0(nkpt0),nband(nkpt*nsppol),nban_dp_rd(nkpt0*nsppol0)
1154  integer, intent(in) :: ngfft(18),npwarr(nkpt),npwarr0(nkpt0),symrel(3,3,nsym)
1155  real(dp), intent(in) :: gmet(3,3),gmet0(3,3),kptns(3,nkpt),kptns0(3,nkpt0),rprimd(3,3)
1156  real(dp), intent(in) :: tnons(3,nsym)
1157  real(dp), intent(out) :: eigen((2*mband)**formeig*mband*nkpt*nsppol)
1158  real(dp), intent(inout) :: cg(2,mcg),cg_disk(2,mcg_disk) !vz_i pw_ortho
1159  real(dp), intent(inout) :: occ(mband*nkpt*nsppol)
1160 
1161 !Local variables-------------------------------
1162  integer :: band_index,band_index_trial,ceksp,debug,dim_eig_k,iband,icg
1163  integer :: icg_disk,icg_trial,idum,ierr,ii,ikassoc,ikassoc_trial,ikpt,ikpt0
1164  integer :: ikpt10,ikpt_trial,ikptsp,ikptsp_old,inplace,isp,isp_max,isppol,isppol0
1165 ! integer :: ipw ! commented out below
1166  integer :: isppol_trial,me,mgfft,my_nspinor,my_nspinor0
1167  integer :: nban_dp_k,nban_dp_rdk,nband_k,nband_rdk,nband_trial,nbd,nbd_max
1168  integer :: ncopy,nkpt_eff,nproc_max,npw0_k,npw_k,npw_ktrial
1169  integer :: read_cg,read_cg_disk,sender,spaceComm
1170  character(len=500) :: message
1171  integer,allocatable :: band_index_k(:,:),icg_k(:,:),kg0_k(:,:),kg_k(:,:)
1172  real(dp) :: tsec(2)
1173  real(dp),allocatable :: eig0_k(:),eig_k(:),occ0_k(:),occ_k(:)
1174 #if defined HAVE_MPI
1175  integer :: iproc,my_ikpt
1176  integer :: tag,test_cycle
1177  integer :: statux(MPI_STATUS_SIZE)
1178  integer,allocatable :: ikassoc_me(:),ikpt_me(:),isppol_me(:),nband_k_me(:)
1179 #endif
1180  integer :: nkpt_max=50
1181 
1182 ! *************************************************************************
1183 
1184 !DEBUG
1185 !write(std_out,*)' wfsinp : enter'
1186 !write(std_out,*)' wfsinp : nband=',nband(:)
1187 !write(std_out,*)' wfsinp : nban_dp_rd=',nban_dp_rd(:)
1188 !write(std_out,*)' wfsinp : localrdwf=',localrdwf
1189 !write(std_out,*)' wfsinp : paralbd,formeig=',mpi_enreg%paralbd,formeig
1190 !write(std_out,*)' wfsinp : indkk0(:,1)=',indkk0(:,1)
1191 !ENDDEBUG
1192 
1193  call timab(720,1,tsec)
1194  call timab(721,3,tsec)
1195 
1196  nkpt_max=50; if(xmpi_paral==1)nkpt_max=-1
1197  nbd_max=size(mpi_enreg%proc_distrb,2)
1198  isp_max=size(mpi_enreg%proc_distrb,3)
1199 
1200 !Init mpi_comm
1201  spaceComm=mpi_enreg%comm_cell
1202  nproc_max=xmpi_comm_size(spaceComm)
1203  me=mpi_enreg%me_kpt
1204  sender = 0
1205 
1206 #if defined HAVE_MPI
1207  if(localrdwf==0)then
1208    ABI_ALLOCATE(ikpt_me,(nproc_max))
1209    ABI_ALLOCATE(nband_k_me,(nproc_max))
1210    ABI_ALLOCATE(ikassoc_me,(nproc_max))
1211    ABI_ALLOCATE(isppol_me,(nproc_max))
1212  end if
1213 #endif
1214 
1215 !Check the validity of formeig
1216  if(formeig/=0.and.formeig/=1)then
1217    write(message, '(a,i0,a)' )' formeig=',formeig,' , but the only allowed values are 0 or 1.'
1218    MSG_BUG(message)
1219  end if
1220 
1221  my_nspinor =max(1,nspinor /mpi_enreg%nproc_spinor)
1222  my_nspinor0=max(1,nspinor0/mpi_enreg%nproc_spinor)
1223 
1224  nkpt_eff=max(nkpt0,nkpt)
1225  if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max
1226 
1227  ABI_ALLOCATE(icg_k,(nkpt,nsppol))
1228  ABI_ALLOCATE(band_index_k,(nkpt,nsppol))
1229 
1230 !write(std_out,*)' wfsinp : me,isppol,ikpt,icg_k,band_index_k'
1231 
1232 !Compute the locations of the blocks in cg, eig and occ
1233  icg=0
1234  band_index=0
1235  ikpt10=0
1236 
1237  do isppol=1,nsppol
1238    do ikpt=1,nkpt
1239 
1240      nband_k=nband(ikpt+(isppol-1)*nkpt)
1241      band_index_k(ikpt,isppol)=band_index
1242 
1243 #if defined HAVE_MPI
1244      test_cycle=0;nbd=min(nband_k,nbd_max);isp=min(isppol,isp_max)
1245      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nbd,isp,me))test_cycle=1
1246      if(test_cycle==1)then
1247        band_index=band_index+nband_k*(2*nband_k)**formeig
1248 !      In the case this k point does not belong to me, cycle
1249        cycle
1250      end if
1251 #endif
1252 
1253      npw_k=npwarr(ikpt)
1254      icg_k(ikpt,isppol)=icg
1255      icg=icg+npw_k*my_nspinor*nband_k
1256 
1257      band_index=band_index+nband_k*(2*nband_k)**formeig
1258 !    write(std_out,'(5i8)' )me,isppol,ikpt,icg_k(ikpt,isppol),band_index_k(ikpt,isppol)
1259    end do ! End k point loop
1260  end do  ! End spin loop
1261 
1262  band_index=0
1263  ikptsp_old=0
1264 
1265 !DEBUG
1266 !write(std_out,*)' wfsinp: before loop'
1267 !write(std_out,*)' nsppol0,nsppol,nkpt0',nsppol0,nsppol,nkpt0
1268 !write(std_out,*)' mpw,mgfft,mpw,mpw0',mpw,mgfft,mpw,mpw0
1269 !ENDDEBUG
1270 
1271  mgfft=maxval(ngfft(1:3))
1272  if(squeeze==1)then
1273    ABI_ALLOCATE(kg_k,(3,mpw))
1274    ABI_ALLOCATE(kg0_k,(3,mpw0))
1275  end if
1276 
1277  eigen(:)=0.0_dp
1278 !occ(:)=0.0_dp
1279 
1280  call timab(721,2,tsec)
1281 
1282 !Loop over spins
1283 !For the time being, do not allow nsppol=2 to nspinor=2 conversion
1284 !MT 20110707: this can be done by a fake call to the routine: see inwffil
1285  do isppol0=1,min(nsppol0,nsppol)
1286 
1287 !  Loop on k points  : get the cg then eventually write on unwfnow
1288    do ikpt0=1,nkpt0
1289 
1290      call timab(722,1,tsec)
1291 
1292      nban_dp_rdk=nban_dp_rd(ikpt0+(isppol0-1)*nkpt0)
1293 
1294 !    DEBUG
1295 !    write(std_out,*)' wfsinp: ikpt0,isppol0,nkpt0=',ikpt0,isppol0,nkpt0
1296 !    write(std_out,*)' nban_dp_rdk=',nban_dp_rdk
1297 !    ENDDEBUG
1298 
1299      npw0_k=npwarr0(ikpt0)
1300      if(ikpt0<=nkpt_eff)then
1301        write(message,'(a,a,2i4)')ch10,' wfsinp: inside loop, init ikpt0,isppol0=',ikpt0,isppol0
1302        call wrtout(std_out,message,'PERS')
1303      end if
1304 
1305 !    Must know whether this k point is needed, and in which
1306 !    block (ikpt, isppol), the wavefunction is to be placed.
1307 !    Select the one for which the number of bands is the biggest.
1308      ikpt=0
1309      isppol=0
1310      ikassoc=0
1311      nband_k=0
1312 #if defined HAVE_MPI
1313      if(localrdwf==0)then
1314        nband_k_me(:)=0
1315        ikpt_me(:)=0
1316        isppol_me(:)=0
1317        ikassoc_me(:)=0
1318        nband_k_me(:)=0
1319      end if
1320 #endif
1321 
1322      do isppol_trial=1,nsppol
1323 
1324        if(nsppol==2 .and. nsppol0==2 .and. isppol0/=isppol_trial)cycle
1325 
1326        do ikassoc_trial=1,nkassoc
1327 
1328          ikpt_trial=indkk0(ikpt0,ikassoc_trial)
1329          if(sppoldbl==2)then
1330            if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle
1331            if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle
1332            if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt
1333          end if
1334 
1335 #if defined HAVE_MPI
1336          if(localrdwf==1)then
1337            if(ikpt_trial/=0)then
1338              nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1339              nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
1340              if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me))ikpt_trial=0
1341            end if
1342          end if
1343 #endif
1344 
1345          if(ikpt_trial/=0)then
1346            nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1347            if(nband_k<nband_trial)then
1348              nband_k=nband_trial ; ikpt=ikpt_trial ; isppol=isppol_trial
1349              ikassoc=ikassoc_trial
1350            end if
1351 
1352 #if defined HAVE_MPI
1353            if(localrdwf==0)then
1354              do iproc=1,nproc_max
1355                my_ikpt=1
1356                nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1357                nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
1358                if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,(iproc-1))) my_ikpt=0
1359                if(my_ikpt/=0)then
1360                  nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1361                  if(nband_k_me(iproc)<nband_trial)then
1362                    nband_k_me(iproc)=nband_trial
1363                    ikpt_me(iproc)=ikpt_trial
1364                    isppol_me(iproc)=isppol_trial
1365                    ikassoc_me(iproc)=ikassoc_trial
1366                  end if
1367                end if
1368              end do
1369            end if
1370 #endif
1371 
1372          end if
1373 
1374        end do ! ikassoc_trial
1375      end do ! isppol_trial
1376 
1377 !    DEBUG
1378 !    write(std_out,*)' wfsinp : me,select isppol,ikpt=',me,isppol,ikpt
1379 #if defined HAVE_MPI
1380 !    write(std_out,*)' wfsinp : me,ikpt_me(:)=',me,ikpt_me(:)
1381 #endif
1382 !    write(std_out,*)' wfsinp : me,isppol_me(:)=',me,isppol_me(:)
1383 !    stop
1384 !    ENDDEBUG
1385 
1386      call timab(722,2,tsec)
1387 
1388 !    If the wavefunction block to be read is interesting ...
1389      if (ikpt/=0)then
1390 
1391        call timab(723,3,tsec)
1392        sender = me
1393        nband_k=nband(ikpt+(isppol-1)*nkpt)
1394        npw_k=npwarr(ikpt)
1395 
1396 #if defined HAVE_MPI
1397        if (localrdwf==1.or.(localrdwf==0.and.me==0)) then
1398 
1399          if(ikpt<=nkpt_eff)then
1400            write(message,'(a,i6,a,i8,a,i4,a,i4)') &
1401 &           ' wfsinp: treating ',nband_k,' bands with npw=',npw_k,' for ikpt=',ikpt,' by node ',me
1402            call wrtout(std_out,message,'PERS')
1403          else if(ikpt==nkpt_eff+1)then
1404            call wrtout(std_out,' wfsinp: prtvol=0 or 1, do not print more k-points.','PERS')
1405          end if
1406 
1407        end if
1408 #endif
1409 
1410        nband_rdk=nban_dp_rdk
1411        if(squeeze==1)nband_rdk=(nban_dp_rdk/nspinor0)*nspinor
1412        if(formeig==0)then
1413          ABI_ALLOCATE(eig_k,(nband_rdk))
1414          ABI_ALLOCATE(occ_k,(nband_rdk))
1415          ABI_ALLOCATE(eig0_k,(nban_dp_rdk))
1416          ABI_ALLOCATE(occ0_k,(nban_dp_rdk))
1417          dim_eig_k=nband_rdk
1418        else if(formeig==1)then
1419          ABI_ALLOCATE(eig_k,(2*nband_rdk*nband_rdk))
1420          ABI_ALLOCATE(eig0_k,(2*nban_dp_rdk*nban_dp_rdk))
1421          dim_eig_k=2*nband_rdk*nband_rdk
1422          ABI_ALLOCATE(occ0_k,(0))
1423          ABI_ALLOCATE(occ_k,(0))
1424        else
1425          ABI_ALLOCATE(occ0_k,(0))
1426          ABI_ALLOCATE(eig0_k,(0))
1427          ABI_ALLOCATE(occ_k,(0))
1428          ABI_ALLOCATE(eig_k,(0))
1429        end if
1430        eig_k(:)=0.0_dp
1431        eig0_k(:)=0.0_dp
1432 
1433 !      Generate or read the cg for this k point
1434 !      Either read into cg, or read into cg_disk
1435        read_cg=1 ; read_cg_disk=0
1436        if(squeeze==1)then
1437          read_cg=0 ; read_cg_disk=1
1438        end if
1439 #if defined HAVE_MPI
1440        if(localrdwf==0)then
1441          read_cg=0
1442          read_cg_disk=0
1443 !        XG20040106 The following condition is correct
1444          if(me==0)read_cg_disk=1
1445        end if
1446 #endif
1447 
1448        icg=0
1449        if(read_cg==1)icg=icg_k(ikpt,isppol)
1450 
1451 !      DEBUG
1452 !      write(std_out,*)' wfsinp: before initwf',wff1%offwff
1453 !      write(std_out,*)' wfsinp: me,read_cg,read_cg_disk=',me,read_cg,read_cg_disk
1454 !      write(std_out,*)' wfsinp: nban_dp_rdk=',nban_dp_rdk
1455 !      ENDDEBUG
1456        call timab(723,2,tsec)
1457        call timab(724,3,tsec)
1458 
1459        if(read_cg_disk==1)then
1460          call initwf (cg_disk,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,&
1461 &         isppol0,mcg_disk,mpi_enreg0, &
1462 &         nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1)
1463        end if
1464 
1465        if(read_cg==1)then
1466          call initwf (cg,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,&
1467 &         isppol0,mcg,mpi_enreg0,&
1468 &         nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1)
1469        end if
1470 
1471        call timab(724,2,tsec)
1472        call timab(725,3,tsec)
1473 
1474        nban_dp_k=min(nban_dp_rdk,(nband_k/nspinor)*nspinor0)
1475 !      This band_index is defined BEFORE the eventual redefinition
1476 !      of ikpt and isppol, needed  when localrdwf==0 in parallel
1477        band_index=band_index_k(ikpt,isppol)
1478 !      DEBUG
1479 !      write(std_out,*)' wfsinp: me,cg_disk(:,1)=',me,cg_disk(:,1)
1480 !      ENDDEBUG
1481 
1482 !      DEBUG
1483 !      if(me==0 .and. ikpt0==1)then
1484 !      write(std_out,*)' wfsinp : cg array, before trial, ikpt0=',ikpt0
1485 !      do ipw=1,15
1486 !      write(std_out,'(i4,2es20.10)' )ipw,cg(:,ipw)
1487 !      end do
1488 !      end if
1489 !      ENDDEBUG
1490 
1491 
1492 #if defined HAVE_MPI
1493        if(localrdwf==0)then
1494 !        Warning: In that case , not yet // on nspinors
1495 !        Transmit to each of the other processors, when needed
1496          if(nproc_max>=2)then
1497            do iproc=2,nproc_max
1498 !            Only me=0 and me=iproc-1 are concerned by this
1499              if(me==0 .or. me==iproc-1)then
1500 
1501                ikpt=ikpt_me(iproc)
1502                isppol=isppol_me(iproc)
1503 
1504                if(ikpt/=0)then
1505 !                In this case, processor iproc-1 needs the data
1506 !                Generate a common tag
1507                  tag=256*(ikpt-1)+iproc+1
1508                  if(isppol==2)tag=-tag
1509                  nband_k=nband(ikpt+(isppol-1)*nkpt)
1510                  npw_k=npwarr(ikpt)
1511 !                SEND
1512                  if(me==0)then
1513                    write(std_out,*)'SENDWFSINP ',me
1514                    call MPI_SEND(cg_disk,2*npw_k*my_nspinor*nband_k,&
1515 &                   MPI_DOUBLE_PRECISION,iproc-1,tag,spaceComm,ierr)
1516                  end if
1517 !                RECEIVE
1518                  if(me==iproc-1)then
1519                    call MPI_RECV(cg_disk,2*npw_k*my_nspinor*nband_k,&
1520 &                   MPI_DOUBLE_PRECISION,0,tag,spaceComm,statux,ierr)
1521                    icg=icg_k(ikpt,isppol)
1522                    if(squeeze==0)then
1523                      cg(:,icg+1:icg+npw_k*my_nspinor*nband_k)=&
1524 &                     cg_disk(:,1:npw_k*my_nspinor*nband_k)
1525                    end if
1526                    ikassoc=ikassoc_me(iproc)
1527                  end if
1528                end if
1529 
1530              end if
1531            end do ! iproc
1532          end if
1533 
1534 !        DEBUG
1535 !        write(std_out,*)' wfsinp: me, iproc loop finished',me
1536 !        ENDDEBUG
1537 
1538 !        Take care of me=0 needing the data
1539          if (me==0) then
1540            ikpt=ikpt_me(me+1)
1541            isppol=isppol_me(me+1)
1542            if(ikpt/=0 )then
1543              nband_k=nband(ikpt+(isppol-1)*nkpt)
1544              npw_k=npwarr(ikpt)
1545 !            I am the master node, and I might need my own data
1546              icg=icg_k(ikpt,isppol)
1547              if(squeeze==0)then
1548 !              Copy from cg_disk to cg
1549                cg(:,1+icg:npw_k*my_nspinor*nband_k+icg)= &
1550 &               cg_disk(:,1:npw_k*my_nspinor*nband_k)
1551              end if
1552              ikassoc=ikassoc_me(me+1)
1553            end if
1554          end if
1555 !        For the eigenvalues and occ, the transmission is much easier to write !
1556          call MPI_BCAST(eig0_k,nban_dp_rdk*(2*nban_dp_rdk)**formeig ,&
1557 &         MPI_DOUBLE_PRECISION,0,spaceComm,ierr)
1558        end if
1559 #endif
1560 
1561        if(formeig==0)then
1562 !        The transfer from eig0_k to eig_k uses nban_dp_rdk, which contains
1563 !        the maximal information.
1564          if(nspinor0==nspinor .or. squeeze==0)then
1565            eig_k(1:nban_dp_rdk)=eig0_k(1:nban_dp_rdk)
1566            occ_k(1:nban_dp_rdk)=occ0_k(1:nban_dp_rdk)
1567          else if(nspinor0==1 .and. nspinor==2)then
1568            do iband=1,nban_dp_rdk
1569              eig_k(2*iband  )=eig0_k(iband)
1570              eig_k(2*iband-1)=eig0_k(iband)
1571              occ_k(2*iband  )=occ0_k(iband)*0.5_dp
1572              occ_k(2*iband-1)=occ0_k(iband)*0.5_dp
1573            end do
1574          else if(nspinor0==2 .and. nspinor==1)then
1575            do iband=1,nban_dp_rdk
1576              eig_k(iband)=eig0_k(2*iband-1)
1577              occ_k(iband)=occ0_k(2*iband-1)*2.0_dp
1578            end do
1579          end if
1580 
1581 !        DEBUG
1582 !        write(std_out,*)' wfsinp: me,band_index,ikpt,isppol',me,band_index,ikpt,isppol
1583 !        ENDDEBUG
1584 
1585 !        The transfer to eigen uses nban_dp_k, that is bound by the number
1586 !        of bands for this k point.
1587          ncopy=min(dim_eig_k,(nban_dp_k/nspinor0)*nspinor)
1588          eigen(1+band_index:ncopy+band_index)=eig_k(1:ncopy)
1589 !        The transfer of occ should be done.
1590 
1591 #if defined HAVE_MPI
1592          if(localrdwf==0 .and. ikpt/=0)then
1593 !          Do not forget : ikpt,isppol were redefined ...
1594            band_index=band_index_k(ikpt,isppol)
1595            eigen(1+band_index:(nban_dp_k/nspinor0)*nspinor+band_index) = eig_k(1:(nban_dp_k/nspinor0)*nspinor)
1596 !          The transfer of occ should be done.
1597          end if
1598 #endif
1599 
1600        else if(formeig==1)then
1601          call wrtout(std_out,ABI_FUNC//': transfer of first-order eigs not yet coded!',"COLL")
1602        end if
1603 
1604 !      DEBUG
1605 !      write(std_out,*)' wfsinp : me,transferred eig_k',me
1606 !      write(std_out,*)' me,mkmem,nsppol,nsppol0,isppol',me,mkmem,nsppol,nsppol0,isppol
1607 !      write(std_out,*)' me,nkassoc,ikassoc',me,nkassoc,ikassoc
1608 !      ENDDEBUG
1609 
1610        call timab(725,2,tsec)
1611 
1612 !      Write to disk if appropriate
1613 !      The coding has to be done ... here, only fragments ...
1614        call timab(727,3,tsec)
1615 
1616        do isppol_trial=1,nsppol
1617          if(nsppol==2 .and. nsppol0==2 .and. isppol_trial/=isppol)cycle
1618          do ikassoc_trial=1,nkassoc
1619 
1620 !            DEBUG
1621 !            write(std_out,*)' wfsinp: me, for ikassoc,isppol',&
1622 !            &       me,ikassoc,isppol
1623 !            write(std_out,*)' wfsinp: me, try ikassoc_trial,isppol_trial,nband_k',&
1624 !            &       me,ikassoc_trial,isppol_trial,nband_k
1625 !            ENDDEBUG
1626 
1627 !            No conversion is to be done : it will be converted in newkpt
1628 !            If squeeze==0, the block with the ikpt corresponding to ikassoc,
1629 !            and with isppol, contains the wavefunction already
1630            if( ikassoc_trial/=ikassoc                   .or. &
1631 &           (isppol_trial/=isppol .and. sppoldbl==1 ).or. &
1632 &           squeeze==1                                       )then
1633 
1634              ikpt_trial=indkk0(ikpt0,ikassoc_trial)
1635              if(sppoldbl==2)then
1636                if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle
1637                if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle
1638                if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt
1639              end if
1640 
1641 
1642 #if defined HAVE_MPI
1643              if(ikpt_trial/=0)then
1644                nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1645                nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
1646                if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me)) ikpt_trial=0
1647              end if
1648 #endif
1649 
1650              if(ikpt_trial/=0 .and. ikpt_trial<=nkpt_eff)then
1651                write(message,'(2a,2i5)')ch10,' wfsinp: transfer to ikpt_trial,isppol_trial=',ikpt_trial,isppol_trial
1652                call wrtout(std_out,message,'PERS')
1653              end if
1654 
1655              if(ikpt_trial/=0)then
1656                icg_trial=icg_k(ikpt_trial,isppol_trial)
1657                band_index_trial=band_index_k(ikpt_trial,isppol_trial)
1658                nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1659                nban_dp_k=min(nban_dp_rdk,(nband_trial/nspinor)*nspinor0)
1660 
1661                if(squeeze==0)then
1662 !                  GMR: modified to avoid compiler bug
1663 !                  cg(:,1+icg_trial:npw0_k*my_nspinor0*nband_trial+icg_trial)=&
1664 !                  &          cg(:,1+icg:npw0_k*my_nspinor0*nband_trial+icg)
1665                  do ii=1,npw0_k*my_nspinor0*nband_trial
1666                    cg(:,ii+icg_trial)=cg(:,ii+icg)
1667                  end do
1668 !                  GMR
1669                  if(formeig==0)then
1670                    eigen(1+band_index_trial:nban_dp_k+band_index_trial)=eig_k(1:nban_dp_k)
1671 !                    occ(1+band_index_trial:nban_dp_k+band_index_trial)=&
1672 !                    &           occ_k(1:nban_dp_k)
1673                  end if
1674 !                  RF transfer of eigenvalues still to be coded
1675                else if(squeeze==1)then
1676                  npw_ktrial=npwarr(ikpt_trial)
1677                  nband_k=(nban_dp_k/nspinor0)*nspinor
1678 !                  Conversion to be done
1679                  ceksp=0 ; debug=0 ; icg_disk=0 ; idum=0 ; inplace=0
1680 !                  Note that this routine also convert eig and occ
1681 !                  even if the conversion had already been done
1682 
1683 
1684                  call wfconv(ceksp,cg_disk,cg,debug,ecut0,ecut,ecut_eff,&
1685 &                 eig0_k,eig_k,exchn2n3d,formeig,gmet0,gmet,&
1686 &                 icg_disk,icg_trial,ikpt0,ikpt10,ikpt_trial,indkk,&
1687 &                 inplace,isppol_trial,istwfk0,istwfk,&
1688 &                 kg0_k,kg_k,kptns0,kptns,nban_dp_rdk,nband_rdk,&
1689 &                 mcg_disk,mcg,mpi_enreg0,mpi_enreg,mpw0,mpw,&
1690 &                 nban_dp_rdk,nband_trial,ngfft,ngfft,nkpt0,nkpt,&
1691 &                 npw0_k,npw_ktrial,nspinor0,nspinor,nsym,&
1692 &                 occ0_k,occ_k,optorth,randalg,restart,rprimd,&
1693 &                 sppoldbl,symrel,tnons)
1694 
1695 !                  DEBUG
1696 !                  write(std_out,*)' wfsinp: ikpt_trial=',ikpt_trial
1697 !                  write(std_out,*)' nband_k,band_index_trial',nband_k,band_index_trial
1698 !                  write(std_out,*)eig0_k(1:nban_dp_k)
1699 !                  write(std_out,*)eig_k(1:nband_k)
1700 !                  ENDDEBUG
1701                  eigen(1+band_index_trial:nband_k+band_index_trial)=eig_k(1:nband_k)
1702 !                  occ(1+band_index_trial:nband_k+band_index_trial)=&
1703 !                  &           occ_k(1:nband_k)
1704 !                  RF transfer of eigenvalues still to be coded
1705 
1706 !                  Endif squeeze==1
1707                end if
1708 
1709 !                DEBUG
1710 !                if(ikpt_trial==2)then
1711 !                write(std_out,*)' wfsinp: iband,ipw,cg for ikpt_trial=2'
1712 !                write(std_out,*)' nband_trial,npw_ktrial=',nband_trial,npw_ktrial
1713 !                do iband=1,nband_trial
1714 !                do ipw=1,npw_ktrial
1715 !                write(std_out,'(2i5,2es16.6)' )&
1716 !                &             iband,ipw,cg(:,ipw+(iband-1)*npw_ktrial+icg_trial)
1717 !                end do
1718 !                end do
1719 !                end if
1720 !                ENDDEBUG
1721 
1722 !                End if ikpt_trial/=0
1723              end if
1724 
1725 !              End if ikpt_trial already initialized
1726            end if
1727 
1728          end do ! ikassoc_trial
1729        end do ! isppol_trial
1730 
1731        call timab(727,2,tsec)
1732 
1733        ABI_DEALLOCATE(eig_k)
1734        ABI_DEALLOCATE(eig0_k)
1735        !if(formeig==0) then
1736        ABI_DEALLOCATE(occ_k)
1737        ABI_DEALLOCATE(occ0_k)
1738        !end if
1739 
1740      end if  ! End the condition of need of this k point
1741    end do ! End of the k loop
1742  end do !  End of spin loop
1743 
1744 #if defined HAVE_MPI
1745  call timab(67,1,tsec)
1746 
1747 !Still need to skip last k points to read history (MS)
1748 !WARNING : not yet for formeig=1, but is it needed ?
1749  if(formeig==0)then
1750    do ikptsp=ikptsp_old+1,nkpt0*nsppol0
1751      isppol=1 ; if(ikptsp>nkpt0)isppol=2
1752      ikpt=ikptsp-nkpt0*(isppol-1)
1753      call WffReadSkipK(formeig,headform0,ikpt,isppol,mpi_enreg,wff1)
1754    end do
1755  end if
1756 
1757 !Transmit eigenvalues. This routine works in both localrdwf=0 or 1 cases.
1758  call pareigocc(eigen,formeig,localrdwf,mpi_enreg,mband,nband,nkpt,nsppol,occ,1)
1759 
1760 
1761  if(localrdwf==0)then
1762    ABI_DEALLOCATE(ikpt_me)
1763    ABI_DEALLOCATE(nband_k_me)
1764    ABI_DEALLOCATE(ikassoc_me)
1765    ABI_DEALLOCATE(isppol_me)
1766  end if
1767 
1768  call timab(67,2,tsec)
1769 #endif
1770 
1771 !****************************************************************************
1772 
1773  if(squeeze==1)then
1774    ABI_DEALLOCATE(kg_k)
1775    ABI_DEALLOCATE(kg0_k)
1776  end if
1777 
1778 
1779 !DEBUG
1780 !if(me==0)then
1781 !write(std_out,*)' wfsinp : cg array='
1782 !icg=0
1783 !do isppol=1,nsppol
1784 !do ikpt=1,1
1785 !nband_k=nband(ikpt+(isppol-1)*nkpt)
1786 !npw_k=npwarr(ikpt)
1787 !do iband=1,nband_k
1788 !write(std_out,*)' new band, icg=',icg
1789 !do ipw=1,npw_k
1790 !write(std_out,'(4i4,2es20.10)' )isppol,ikpt,iband,ipw,cg(:,icg+ipw)
1791 !end do
1792 !icg=icg+npw_k
1793 !end do
1794 !end do
1795 !end do
1796 !end if
1797 !if(ireadwf==1)stop
1798 !write(std_out,*)' wfsinp : eigen array='
1799 !do ikpt=1,nkpt
1800 !do iband=1,mband
1801 !write(std_out,*)'ikpt,iband,eigen',ikpt,iband,eigen(iband+(ikpt-1)*mband)
1802 !end do
1803 !end do
1804 !ENDDEBUG
1805 
1806  ABI_DEALLOCATE(icg_k)
1807  ABI_DEALLOCATE(band_index_k)
1808 
1809  call timab(720,2,tsec)
1810 
1811 end subroutine wfsinp