TABLE OF CONTENTS
ABINIT/wfsinp [ 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)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, AR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
ecut0=kinetic energy cutoffs for basis sphere 0 (hartree) (if squeeze=1) ecut=kinetic energy cutoffs beyond which the coefficients of cg vanish (Ha) (needed only if squeeze=1) ecut_eff=effective kinetic energy planewave cutoff (hartree), needed to generate the sphere of plane wave 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=informations about MPI parallelization mpi_enreg0=informations 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 informations for input and output files
OUTPUT
if ground state format (formeig=0):
eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
if respfn format (formeig=1):
eigen(2*mband*mband*nkpt*nsppol)=
matrix of eigenvalues (input or init to large number), (Ha)
Conditional output:
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
116 #if defined HAVE_CONFIG_H 117 #include "config.h" 118 #endif 119 120 #include "abi_common.h" 121 122 123 subroutine wfsinp(cg,cg_disk,ecut,ecut0,ecut_eff,eigen,exchn2n3d,& 124 & formeig,gmet,gmet0,headform0,indkk,indkk0,istwfk,& 125 & istwfk0,kptns,kptns0,localrdwf,mband,& 126 & mcg,mcg_disk,mpi_enreg,mpi_enreg0,mpw,mpw0,nband,nban_dp_rd,& 127 & ngfft,nkassoc,nkpt,nkpt0,npwarr,npwarr0,nspinor,& 128 & nspinor0,nsppol,nsppol0,nsym,occ,optorth,prtvol,randalg,restart,rprimd,& 129 & sppoldbl,squeeze,symrel,tnons,wff1) 130 131 use defs_basis 132 use defs_abitypes 133 use m_errors 134 use m_profiling_abi 135 use m_wffile 136 use m_xmpi 137 #if defined HAVE_MPI2 138 use mpi 139 #endif 140 141 use m_occ, only : pareigocc 142 143 !This section has been created automatically by the script Abilint (TD). 144 !Do not modify the following lines by hand. 145 #undef ABI_FUNC 146 #define ABI_FUNC 'wfsinp' 147 use interfaces_14_hidewrite 148 use interfaces_18_timing 149 use interfaces_32_util 150 use interfaces_62_iowfdenpot 151 use interfaces_66_wfs 152 !End of the abilint section 153 154 implicit none 155 156 #if defined HAVE_MPI1 157 include 'mpif.h' 158 #endif 159 160 !Arguments ------------------------------------ 161 integer, intent(in) :: exchn2n3d,formeig,headform0,localrdwf,mband,mcg,mcg_disk 162 integer, intent(in) :: mpw,mpw0,nkassoc,nkpt,nkpt0,nspinor,nspinor0,nsppol,nsppol0,nsym 163 integer, intent(in) :: optorth,prtvol,randalg,restart,sppoldbl,squeeze 164 real(dp), intent(in) :: ecut,ecut0,ecut_eff 165 type(MPI_type), intent(inout) :: mpi_enreg,mpi_enreg0 166 type(wffile_type), intent(inout) :: wff1 167 integer, intent(in) :: indkk(nkpt*sppoldbl,6),indkk0(nkpt0,nkassoc),istwfk(nkpt) 168 integer, intent(in) :: istwfk0(nkpt0),nband(nkpt*nsppol),nban_dp_rd(nkpt0*nsppol0) 169 integer, intent(in) :: ngfft(18),npwarr(nkpt),npwarr0(nkpt0),symrel(3,3,nsym) 170 real(dp), intent(in) :: gmet(3,3),gmet0(3,3),kptns(3,nkpt),kptns0(3,nkpt0),rprimd(3,3) 171 real(dp), intent(in) :: tnons(3,nsym) 172 real(dp), intent(out) :: eigen((2*mband)**formeig*mband*nkpt*nsppol) 173 real(dp), intent(inout) :: cg(2,mcg),cg_disk(2,mcg_disk) !vz_i pw_ortho 174 real(dp), intent(inout) :: occ(mband*nkpt*nsppol) 175 176 !Local variables------------------------------- 177 integer :: band_index,band_index_trial,ceksp,debug,dim_eig_k,iband,icg 178 integer :: icg_disk,icg_trial,idum,ierr,ii,ikassoc,ikassoc_trial,ikpt,ikpt0 179 integer :: ikpt10,ikpt_trial,ikptsp,ikptsp_old,inplace,isp,isp_max,isppol,isppol0 180 ! integer :: ipw ! commented out below 181 integer :: isppol_trial,me,mgfft,my_nspinor,my_nspinor0 182 integer :: nban_dp_k,nban_dp_rdk,nband_k,nband_rdk,nband_trial,nbd,nbd_max 183 integer :: ncopy,nkpt_eff,nproc_max,npw0_k,npw_k,npw_ktrial 184 integer :: read_cg,read_cg_disk,sender,spaceComm 185 character(len=500) :: message 186 integer,allocatable :: band_index_k(:,:),icg_k(:,:),kg0_k(:,:),kg_k(:,:) 187 real(dp) :: tsec(2) 188 real(dp),allocatable :: eig0_k(:),eig_k(:),occ0_k(:),occ_k(:) 189 #if defined HAVE_MPI 190 integer :: iproc,my_ikpt 191 integer :: tag,test_cycle 192 integer :: statux(MPI_STATUS_SIZE) 193 integer,allocatable :: ikassoc_me(:),ikpt_me(:),isppol_me(:),nband_k_me(:) 194 #endif 195 integer :: nkpt_max=50 196 197 ! ************************************************************************* 198 199 !DEBUG 200 !write(std_out,*)' wfsinp : enter' 201 !write(std_out,*)' wfsinp : nband=',nband(:) 202 !write(std_out,*)' wfsinp : nban_dp_rd=',nban_dp_rd(:) 203 !write(std_out,*)' wfsinp : localrdwf=',localrdwf 204 !write(std_out,*)' wfsinp : paralbd,formeig=',mpi_enreg%paralbd,formeig 205 !write(std_out,*)' wfsinp : indkk0(:,1)=',indkk0(:,1) 206 !ENDDEBUG 207 208 call timab(720,1,tsec) 209 call timab(721,3,tsec) 210 211 nkpt_max=50; if(xmpi_paral==1)nkpt_max=-1 212 nbd_max=size(mpi_enreg%proc_distrb,2) 213 isp_max=size(mpi_enreg%proc_distrb,3) 214 215 !Init mpi_comm 216 spaceComm=mpi_enreg%comm_cell 217 nproc_max=xmpi_comm_size(spaceComm) 218 me=mpi_enreg%me_kpt 219 sender = 0 220 221 #if defined HAVE_MPI 222 if(localrdwf==0)then 223 ABI_ALLOCATE(ikpt_me,(nproc_max)) 224 ABI_ALLOCATE(nband_k_me,(nproc_max)) 225 ABI_ALLOCATE(ikassoc_me,(nproc_max)) 226 ABI_ALLOCATE(isppol_me,(nproc_max)) 227 end if 228 #endif 229 230 !Check the validity of formeig 231 if(formeig/=0.and.formeig/=1)then 232 write(message, '(a,i0,a)' )' formeig=',formeig,' , but the only allowed values are 0 or 1.' 233 MSG_BUG(message) 234 end if 235 236 my_nspinor =max(1,nspinor /mpi_enreg%nproc_spinor) 237 my_nspinor0=max(1,nspinor0/mpi_enreg%nproc_spinor) 238 239 nkpt_eff=max(nkpt0,nkpt) 240 if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max 241 242 ABI_ALLOCATE(icg_k,(nkpt,nsppol)) 243 ABI_ALLOCATE(band_index_k,(nkpt,nsppol)) 244 245 !write(std_out,*)' wfsinp : me,isppol,ikpt,icg_k,band_index_k' 246 247 !Compute the locations of the blocks in cg, eig and occ 248 icg=0 249 band_index=0 250 ikpt10=0 251 252 do isppol=1,nsppol 253 do ikpt=1,nkpt 254 255 nband_k=nband(ikpt+(isppol-1)*nkpt) 256 band_index_k(ikpt,isppol)=band_index 257 258 #if defined HAVE_MPI 259 test_cycle=0;nbd=min(nband_k,nbd_max);isp=min(isppol,isp_max) 260 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nbd,isp,me))test_cycle=1 261 if(test_cycle==1)then 262 band_index=band_index+nband_k*(2*nband_k)**formeig 263 ! In the case this k point does not belong to me, cycle 264 cycle 265 end if 266 #endif 267 268 npw_k=npwarr(ikpt) 269 icg_k(ikpt,isppol)=icg 270 icg=icg+npw_k*my_nspinor*nband_k 271 272 band_index=band_index+nband_k*(2*nband_k)**formeig 273 ! write(std_out,'(5i8)' )me,isppol,ikpt,icg_k(ikpt,isppol),band_index_k(ikpt,isppol) 274 end do ! End k point loop 275 end do ! End spin loop 276 277 band_index=0 278 ikptsp_old=0 279 280 !DEBUG 281 !write(std_out,*)' wfsinp: before loop' 282 !write(std_out,*)' nsppol0,nsppol,nkpt0',nsppol0,nsppol,nkpt0 283 !write(std_out,*)' mpw,mgfft,mpw,mpw0',mpw,mgfft,mpw,mpw0 284 !ENDDEBUG 285 286 mgfft=maxval(ngfft(1:3)) 287 if(squeeze==1)then 288 ABI_ALLOCATE(kg_k,(3,mpw)) 289 ABI_ALLOCATE(kg0_k,(3,mpw0)) 290 end if 291 292 eigen(:)=0.0_dp 293 !occ(:)=0.0_dp 294 295 call timab(721,2,tsec) 296 297 !Loop over spins 298 !For the time being, do not allow nsppol=2 to nspinor=2 conversion 299 !MT 20110707: this can be done by a fake call to the routine: see inwffil 300 do isppol0=1,min(nsppol0,nsppol) 301 302 ! Loop on k points : get the cg then eventually write on unwfnow 303 do ikpt0=1,nkpt0 304 305 call timab(722,1,tsec) 306 307 nban_dp_rdk=nban_dp_rd(ikpt0+(isppol0-1)*nkpt0) 308 309 ! DEBUG 310 ! write(std_out,*)' wfsinp: ikpt0,isppol0,nkpt0=',ikpt0,isppol0,nkpt0 311 ! write(std_out,*)' nban_dp_rdk=',nban_dp_rdk 312 ! ENDDEBUG 313 314 npw0_k=npwarr0(ikpt0) 315 if(ikpt0<=nkpt_eff)then 316 write(message,'(a,a,2i4)')ch10,' wfsinp: inside loop, init ikpt0,isppol0=',ikpt0,isppol0 317 call wrtout(std_out,message,'PERS') 318 end if 319 320 ! Must know whether this k point is needed, and in which 321 ! block (ikpt, isppol), the wavefunction is to be placed. 322 ! Select the one for which the number of bands is the biggest. 323 ikpt=0 324 isppol=0 325 ikassoc=0 326 nband_k=0 327 #if defined HAVE_MPI 328 if(localrdwf==0)then 329 nband_k_me(:)=0 330 ikpt_me(:)=0 331 isppol_me(:)=0 332 ikassoc_me(:)=0 333 nband_k_me(:)=0 334 end if 335 #endif 336 337 do isppol_trial=1,nsppol 338 339 if(nsppol==2 .and. nsppol0==2 .and. isppol0/=isppol_trial)cycle 340 341 do ikassoc_trial=1,nkassoc 342 343 ikpt_trial=indkk0(ikpt0,ikassoc_trial) 344 if(sppoldbl==2)then 345 if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle 346 if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle 347 if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt 348 end if 349 350 #if defined HAVE_MPI 351 if(localrdwf==1)then 352 if(ikpt_trial/=0)then 353 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 354 nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max) 355 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me))ikpt_trial=0 356 end if 357 end if 358 #endif 359 360 if(ikpt_trial/=0)then 361 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 362 if(nband_k<nband_trial)then 363 nband_k=nband_trial ; ikpt=ikpt_trial ; isppol=isppol_trial 364 ikassoc=ikassoc_trial 365 end if 366 367 #if defined HAVE_MPI 368 if(localrdwf==0)then 369 do iproc=1,nproc_max 370 my_ikpt=1 371 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 372 nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max) 373 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,(iproc-1))) my_ikpt=0 374 if(my_ikpt/=0)then 375 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 376 if(nband_k_me(iproc)<nband_trial)then 377 nband_k_me(iproc)=nband_trial 378 ikpt_me(iproc)=ikpt_trial 379 isppol_me(iproc)=isppol_trial 380 ikassoc_me(iproc)=ikassoc_trial 381 end if 382 end if 383 end do 384 end if 385 #endif 386 387 end if 388 389 end do ! ikassoc_trial 390 end do ! isppol_trial 391 392 ! DEBUG 393 ! write(std_out,*)' wfsinp : me,select isppol,ikpt=',me,isppol,ikpt 394 #if defined HAVE_MPI 395 ! write(std_out,*)' wfsinp : me,ikpt_me(:)=',me,ikpt_me(:) 396 #endif 397 ! write(std_out,*)' wfsinp : me,isppol_me(:)=',me,isppol_me(:) 398 ! stop 399 ! ENDDEBUG 400 401 call timab(722,2,tsec) 402 403 ! If the wavefunction block to be read is interesting ... 404 if (ikpt/=0)then 405 406 call timab(723,3,tsec) 407 sender = me 408 nband_k=nband(ikpt+(isppol-1)*nkpt) 409 npw_k=npwarr(ikpt) 410 411 #if defined HAVE_MPI 412 if (localrdwf==1.or.(localrdwf==0.and.me==0)) then 413 414 if(ikpt<=nkpt_eff)then 415 write(message,'(a,i6,a,i8,a,i4,a,i4)') & 416 & ' wfsinp: treating ',nband_k,' bands with npw=',npw_k,' for ikpt=',ikpt,' by node ',me 417 call wrtout(std_out,message,'PERS') 418 else if(ikpt==nkpt_eff+1)then 419 call wrtout(std_out,' wfsinp: prtvol=0 or 1, do not print more k-points.','PERS') 420 end if 421 422 end if 423 #endif 424 425 nband_rdk=nban_dp_rdk 426 if(squeeze==1)nband_rdk=(nban_dp_rdk/nspinor0)*nspinor 427 if(formeig==0)then 428 ABI_ALLOCATE(eig_k,(nband_rdk)) 429 ABI_ALLOCATE(occ_k,(nband_rdk)) 430 ABI_ALLOCATE(eig0_k,(nban_dp_rdk)) 431 ABI_ALLOCATE(occ0_k,(nban_dp_rdk)) 432 dim_eig_k=nband_rdk 433 else if(formeig==1)then 434 ABI_ALLOCATE(eig_k,(2*nband_rdk*nband_rdk)) 435 ABI_ALLOCATE(eig0_k,(2*nban_dp_rdk*nban_dp_rdk)) 436 dim_eig_k=2*nband_rdk*nband_rdk 437 ABI_ALLOCATE(occ0_k,(0)) 438 ABI_ALLOCATE(occ_k,(0)) 439 else 440 ABI_ALLOCATE(occ0_k,(0)) 441 ABI_ALLOCATE(eig0_k,(0)) 442 ABI_ALLOCATE(occ_k,(0)) 443 ABI_ALLOCATE(eig_k,(0)) 444 end if 445 eig_k(:)=0.0_dp 446 eig0_k(:)=0.0_dp 447 448 ! Generate or read the cg for this k point 449 ! Either read into cg, or read into cg_disk 450 read_cg=1 ; read_cg_disk=0 451 if(squeeze==1)then 452 read_cg=0 ; read_cg_disk=1 453 end if 454 #if defined HAVE_MPI 455 if(localrdwf==0)then 456 read_cg=0 457 read_cg_disk=0 458 ! XG20040106 The following condition is correct 459 if(me==0)read_cg_disk=1 460 end if 461 #endif 462 463 icg=0 464 if(read_cg==1)icg=icg_k(ikpt,isppol) 465 466 ! DEBUG 467 ! write(std_out,*)' wfsinp: before initwf',wff1%offwff 468 ! write(std_out,*)' wfsinp: me,read_cg,read_cg_disk=',me,read_cg,read_cg_disk 469 ! write(std_out,*)' wfsinp: nban_dp_rdk=',nban_dp_rdk 470 ! ENDDEBUG 471 call timab(723,2,tsec) 472 call timab(724,3,tsec) 473 474 if(read_cg_disk==1)then 475 call initwf (cg_disk,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,& 476 & isppol0,mcg_disk,mpi_enreg0, & 477 & nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1) 478 end if 479 480 if(read_cg==1)then 481 call initwf (cg,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,& 482 & isppol0,mcg,mpi_enreg0,& 483 & nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1) 484 end if 485 486 call timab(724,2,tsec) 487 call timab(725,3,tsec) 488 489 nban_dp_k=min(nban_dp_rdk,(nband_k/nspinor)*nspinor0) 490 ! This band_index is defined BEFORE the eventual redefinition 491 ! of ikpt and isppol, needed when localrdwf==0 in parallel 492 band_index=band_index_k(ikpt,isppol) 493 ! DEBUG 494 ! write(std_out,*)' wfsinp: me,cg_disk(:,1)=',me,cg_disk(:,1) 495 ! ENDDEBUG 496 497 ! DEBUG 498 ! if(me==0 .and. ikpt0==1)then 499 ! write(std_out,*)' wfsinp : cg array, before trial, ikpt0=',ikpt0 500 ! do ipw=1,15 501 ! write(std_out,'(i4,2es20.10)' )ipw,cg(:,ipw) 502 ! end do 503 ! end if 504 ! ENDDEBUG 505 506 507 #if defined HAVE_MPI 508 if(localrdwf==0)then 509 ! Warning: In that case , not yet // on nspinors 510 ! Transmit to each of the other processors, when needed 511 if(nproc_max>=2)then 512 do iproc=2,nproc_max 513 ! Only me=0 and me=iproc-1 are concerned by this 514 if(me==0 .or. me==iproc-1)then 515 516 ikpt=ikpt_me(iproc) 517 isppol=isppol_me(iproc) 518 519 if(ikpt/=0)then 520 ! In this case, processor iproc-1 needs the data 521 ! Generate a common tag 522 tag=256*(ikpt-1)+iproc+1 523 if(isppol==2)tag=-tag 524 nband_k=nband(ikpt+(isppol-1)*nkpt) 525 npw_k=npwarr(ikpt) 526 ! SEND 527 if(me==0)then 528 write(std_out,*)'SENDWFSINP ',me 529 call MPI_SEND(cg_disk,2*npw_k*my_nspinor*nband_k,& 530 & MPI_DOUBLE_PRECISION,iproc-1,tag,spaceComm,ierr) 531 end if 532 ! RECEIVE 533 if(me==iproc-1)then 534 call MPI_RECV(cg_disk,2*npw_k*my_nspinor*nband_k,& 535 & MPI_DOUBLE_PRECISION,0,tag,spaceComm,statux,ierr) 536 icg=icg_k(ikpt,isppol) 537 if(squeeze==0)then 538 cg(:,icg+1:icg+npw_k*my_nspinor*nband_k)=& 539 & cg_disk(:,1:npw_k*my_nspinor*nband_k) 540 end if 541 ikassoc=ikassoc_me(iproc) 542 end if 543 end if 544 545 end if 546 end do ! iproc 547 end if 548 549 ! DEBUG 550 ! write(std_out,*)' wfsinp: me, iproc loop finished',me 551 ! ENDDEBUG 552 553 ! Take care of me=0 needing the data 554 if (me==0) then 555 ikpt=ikpt_me(me+1) 556 isppol=isppol_me(me+1) 557 if(ikpt/=0 )then 558 nband_k=nband(ikpt+(isppol-1)*nkpt) 559 npw_k=npwarr(ikpt) 560 ! I am the master node, and I might need my own data 561 icg=icg_k(ikpt,isppol) 562 if(squeeze==0)then 563 ! Copy from cg_disk to cg 564 cg(:,1+icg:npw_k*my_nspinor*nband_k+icg)= & 565 & cg_disk(:,1:npw_k*my_nspinor*nband_k) 566 end if 567 ikassoc=ikassoc_me(me+1) 568 end if 569 end if 570 ! For the eigenvalues and occ, the transmission is much easier to write ! 571 call MPI_BCAST(eig0_k,nban_dp_rdk*(2*nban_dp_rdk)**formeig ,& 572 & MPI_DOUBLE_PRECISION,0,spaceComm,ierr) 573 end if 574 #endif 575 576 if(formeig==0)then 577 ! The transfer from eig0_k to eig_k uses nban_dp_rdk, which contains 578 ! the maximal information. 579 if(nspinor0==nspinor .or. squeeze==0)then 580 eig_k(1:nban_dp_rdk)=eig0_k(1:nban_dp_rdk) 581 occ_k(1:nban_dp_rdk)=occ0_k(1:nban_dp_rdk) 582 else if(nspinor0==1 .and. nspinor==2)then 583 do iband=1,nban_dp_rdk 584 eig_k(2*iband )=eig0_k(iband) 585 eig_k(2*iband-1)=eig0_k(iband) 586 occ_k(2*iband )=occ0_k(iband)*0.5_dp 587 occ_k(2*iband-1)=occ0_k(iband)*0.5_dp 588 end do 589 else if(nspinor0==2 .and. nspinor==1)then 590 do iband=1,nban_dp_rdk 591 eig_k(iband)=eig0_k(2*iband-1) 592 occ_k(iband)=occ0_k(2*iband-1)*2.0_dp 593 end do 594 end if 595 596 ! DEBUG 597 ! write(std_out,*)' wfsinp: me,band_index,ikpt,isppol',me,band_index,ikpt,isppol 598 ! ENDDEBUG 599 600 ! The transfer to eigen uses nban_dp_k, that is bound by the number 601 ! of bands for this k point. 602 ncopy=min(dim_eig_k,(nban_dp_k/nspinor0)*nspinor) 603 eigen(1+band_index:ncopy+band_index)=eig_k(1:ncopy) 604 ! The transfer of occ should be done. 605 606 #if defined HAVE_MPI 607 if(localrdwf==0 .and. ikpt/=0)then 608 ! Do not forget : ikpt,isppol were redefined ... 609 band_index=band_index_k(ikpt,isppol) 610 eigen(1+band_index:(nban_dp_k/nspinor0)*nspinor+band_index) = eig_k(1:(nban_dp_k/nspinor0)*nspinor) 611 ! The transfer of occ should be done. 612 end if 613 #endif 614 615 else if(formeig==1)then 616 call wrtout(std_out,ABI_FUNC//': transfer of first-order eigs not yet coded!',"COLL") 617 end if 618 619 ! DEBUG 620 ! write(std_out,*)' wfsinp : me,transferred eig_k',me 621 ! write(std_out,*)' me,mkmem,nsppol,nsppol0,isppol',me,mkmem,nsppol,nsppol0,isppol 622 ! write(std_out,*)' me,nkassoc,ikassoc',me,nkassoc,ikassoc 623 ! ENDDEBUG 624 625 call timab(725,2,tsec) 626 627 ! Write to disk if appropriate 628 ! The coding has to be done ... here, only fragments ... 629 call timab(727,3,tsec) 630 631 do isppol_trial=1,nsppol 632 if(nsppol==2 .and. nsppol0==2 .and. isppol_trial/=isppol)cycle 633 do ikassoc_trial=1,nkassoc 634 635 ! DEBUG 636 ! write(std_out,*)' wfsinp: me, for ikassoc,isppol',& 637 ! & me,ikassoc,isppol 638 ! write(std_out,*)' wfsinp: me, try ikassoc_trial,isppol_trial,nband_k',& 639 ! & me,ikassoc_trial,isppol_trial,nband_k 640 ! ENDDEBUG 641 642 ! No conversion is to be done : it will be converted in newkpt 643 ! If squeeze==0, the block with the ikpt corresponding to ikassoc, 644 ! and with isppol, contains the wavefunction already 645 if( ikassoc_trial/=ikassoc .or. & 646 & (isppol_trial/=isppol .and. sppoldbl==1 ).or. & 647 & squeeze==1 )then 648 649 ikpt_trial=indkk0(ikpt0,ikassoc_trial) 650 if(sppoldbl==2)then 651 if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle 652 if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle 653 if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt 654 end if 655 656 657 #if defined HAVE_MPI 658 if(ikpt_trial/=0)then 659 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 660 nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max) 661 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me)) ikpt_trial=0 662 end if 663 #endif 664 665 if(ikpt_trial/=0 .and. ikpt_trial<=nkpt_eff)then 666 write(message,'(2a,2i5)')ch10,' wfsinp: transfer to ikpt_trial,isppol_trial=',ikpt_trial,isppol_trial 667 call wrtout(std_out,message,'PERS') 668 end if 669 670 if(ikpt_trial/=0)then 671 icg_trial=icg_k(ikpt_trial,isppol_trial) 672 band_index_trial=band_index_k(ikpt_trial,isppol_trial) 673 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 674 nban_dp_k=min(nban_dp_rdk,(nband_trial/nspinor)*nspinor0) 675 676 if(squeeze==0)then 677 ! GMR: modified to avoid compiler bug 678 ! cg(:,1+icg_trial:npw0_k*my_nspinor0*nband_trial+icg_trial)=& 679 ! & cg(:,1+icg:npw0_k*my_nspinor0*nband_trial+icg) 680 do ii=1,npw0_k*my_nspinor0*nband_trial 681 cg(:,ii+icg_trial)=cg(:,ii+icg) 682 end do 683 ! GMR 684 if(formeig==0)then 685 eigen(1+band_index_trial:nban_dp_k+band_index_trial)=eig_k(1:nban_dp_k) 686 ! occ(1+band_index_trial:nban_dp_k+band_index_trial)=& 687 ! & occ_k(1:nban_dp_k) 688 end if 689 ! RF transfer of eigenvalues still to be coded 690 else if(squeeze==1)then 691 npw_ktrial=npwarr(ikpt_trial) 692 nband_k=(nban_dp_k/nspinor0)*nspinor 693 ! Conversion to be done 694 ceksp=0 ; debug=0 ; icg_disk=0 ; idum=0 ; inplace=0 695 ! Note that this routine also convert eig and occ 696 ! even if the conversion had already been done 697 698 699 call wfconv(ceksp,cg_disk,cg,debug,ecut0,ecut,ecut_eff,& 700 & eig0_k,eig_k,exchn2n3d,formeig,gmet0,gmet,& 701 & icg_disk,icg_trial,ikpt0,ikpt10,ikpt_trial,indkk,& 702 & inplace,isppol_trial,istwfk0,istwfk,& 703 & kg0_k,kg_k,kptns0,kptns,nban_dp_rdk,nband_rdk,& 704 & mcg_disk,mcg,mpi_enreg0,mpi_enreg,mpw0,mpw,& 705 & nban_dp_rdk,nband_trial,ngfft,ngfft,nkpt0,nkpt,& 706 & npw0_k,npw_ktrial,nspinor0,nspinor,nsym,& 707 & occ0_k,occ_k,optorth,randalg,restart,rprimd,& 708 & sppoldbl,symrel,tnons) 709 710 ! DEBUG 711 ! write(std_out,*)' wfsinp: ikpt_trial=',ikpt_trial 712 ! write(std_out,*)' nband_k,band_index_trial',nband_k,band_index_trial 713 ! write(std_out,*)eig0_k(1:nban_dp_k) 714 ! write(std_out,*)eig_k(1:nband_k) 715 ! ENDDEBUG 716 eigen(1+band_index_trial:nband_k+band_index_trial)=eig_k(1:nband_k) 717 ! occ(1+band_index_trial:nband_k+band_index_trial)=& 718 ! & occ_k(1:nband_k) 719 ! RF transfer of eigenvalues still to be coded 720 721 ! Endif squeeze==1 722 end if 723 724 ! DEBUG 725 ! if(ikpt_trial==2)then 726 ! write(std_out,*)' wfsinp: iband,ipw,cg for ikpt_trial=2' 727 ! write(std_out,*)' nband_trial,npw_ktrial=',nband_trial,npw_ktrial 728 ! do iband=1,nband_trial 729 ! do ipw=1,npw_ktrial 730 ! write(std_out,'(2i5,2es16.6)' )& 731 ! & iband,ipw,cg(:,ipw+(iband-1)*npw_ktrial+icg_trial) 732 ! end do 733 ! end do 734 ! end if 735 ! ENDDEBUG 736 737 ! End if ikpt_trial/=0 738 end if 739 740 ! End if ikpt_trial already initialized 741 end if 742 743 end do ! ikassoc_trial 744 end do ! isppol_trial 745 746 call timab(727,2,tsec) 747 748 ABI_DEALLOCATE(eig_k) 749 ABI_DEALLOCATE(eig0_k) 750 !if(formeig==0) then 751 ABI_DEALLOCATE(occ_k) 752 ABI_DEALLOCATE(occ0_k) 753 !end if 754 755 end if ! End the condition of need of this k point 756 end do ! End of the k loop 757 end do ! End of spin loop 758 759 #if defined HAVE_MPI 760 call timab(67,1,tsec) 761 762 !Still need to skip last k points to read history (MS) 763 !WARNING : not yet for formeig=1, but is it needed ? 764 if(formeig==0)then 765 do ikptsp=ikptsp_old+1,nkpt0*nsppol0 766 isppol=1 ; if(ikptsp>nkpt0)isppol=2 767 ikpt=ikptsp-nkpt0*(isppol-1) 768 call WffReadSkipK(formeig,headform0,ikpt,isppol,mpi_enreg,wff1) 769 end do 770 end if 771 772 !Transmit eigenvalues. This routine works in both localrdwf=0 or 1 cases. 773 call pareigocc(eigen,formeig,localrdwf,mpi_enreg,mband,nband,nkpt,nsppol,occ,1) 774 775 776 if(localrdwf==0)then 777 ABI_DEALLOCATE(ikpt_me) 778 ABI_DEALLOCATE(nband_k_me) 779 ABI_DEALLOCATE(ikassoc_me) 780 ABI_DEALLOCATE(isppol_me) 781 end if 782 783 call timab(67,2,tsec) 784 #endif 785 786 !**************************************************************************** 787 788 if(squeeze==1)then 789 ABI_DEALLOCATE(kg_k) 790 ABI_DEALLOCATE(kg0_k) 791 end if 792 793 794 !DEBUG 795 !if(me==0)then 796 !write(std_out,*)' wfsinp : cg array=' 797 !icg=0 798 !do isppol=1,nsppol 799 !do ikpt=1,1 800 !nband_k=nband(ikpt+(isppol-1)*nkpt) 801 !npw_k=npwarr(ikpt) 802 !do iband=1,nband_k 803 !write(std_out,*)' new band, icg=',icg 804 !do ipw=1,npw_k 805 !write(std_out,'(4i4,2es20.10)' )isppol,ikpt,iband,ipw,cg(:,icg+ipw) 806 !end do 807 !icg=icg+npw_k 808 !end do 809 !end do 810 !end do 811 !end if 812 !if(ireadwf==1)stop 813 !write(std_out,*)' wfsinp : eigen array=' 814 !do ikpt=1,nkpt 815 !do iband=1,mband 816 !write(std_out,*)'ikpt,iband,eigen',ikpt,iband,eigen(iband+(ikpt-1)*mband) 817 !end do 818 !end do 819 !ENDDEBUG 820 821 ABI_DEALLOCATE(icg_k) 822 ABI_DEALLOCATE(band_index_k) 823 824 call timab(720,2,tsec) 825 826 end subroutine wfsinp