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