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