TABLE OF CONTENTS
ABINIT/m_wfk_analyze [ Modules ]
NAME
m_wfk_analyze
FUNCTION
Post-processing tools for WFK file
COPYRIGHT
Copyright (C) 2008-2022 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_wfk_analyze 23 24 use defs_basis 25 use m_abicore 26 use m_xmpi 27 use m_errors 28 use m_hdr 29 use m_crystal 30 use m_ebands 31 use m_nctk 32 use m_wfk 33 use m_wfd 34 use m_dtset 35 use m_dtfil 36 use m_distribfft 37 38 use defs_datatypes, only : pseudopotential_type, ebands_t 39 use defs_abitypes, only : mpi_type 40 use m_time, only : timab 41 use m_fstrings, only : strcat, sjoin, itoa, ftoa 42 use m_fftcore, only : print_ngfft 43 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 44 use m_esymm, only : esymm_t, esymm_free 45 use m_ddk, only : ddkstore_t 46 use m_pawang, only : pawang_type 47 use m_pawrad, only : pawrad_type 48 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 49 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 50 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 51 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init, pawfgrtab_print 52 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_inquire_dim 53 use m_pawdij, only : pawdij, symdij 54 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 55 use m_paw_sphharm, only : setsym_ylm 56 use m_paw_init, only : pawinit, paw_gencond 57 use m_paw_nhat, only : nhatgrid 58 use m_paw_tools, only : chkpawovlp 59 use m_paw_correlations,only : pawpuxinit 60 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 61 use m_classify_bands, only : classify_bands 62 use m_pspini, only : pspini 63 use m_sigtk, only : sigtk_kpts_in_erange 64 use m_iowf, only : prtkbff 65 66 implicit none 67 68 private
ABINIT/wfk_analyze [ Functions ]
NAME
wfk_analyze
FUNCTION
Main routine implementing postprocessing tools for the WFK file. Main differences wrt cut3d: - MPI support. - No interactive prompt. - Run the analysis once, store all the important results in netcdf files and use python tools to analyze data
INPUTS
acell(3)=Length scales of primitive translations (bohr) codvsn=Code version dtfil<datafiles_type>=Variables related to files. dtset<dataset_type>=All input variables for this dataset. pawang<pawang_type)>=PAW angular mesh and related data. pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data. pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. Before entering the first time in the routine, a significant part of psps has been initialized : the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp. All the remaining components of psps are to be initialized in the call to pspini. The next time the code enters bethe_salpeter, psps might be identical to the one of the previous dtset, in which case, no reinitialisation is scheduled in pspini.F90. rprim(3,3)=Dimensionless real space primitive translations. xred(3,natom)=Reduced atomic coordinates.
NOTES
ON THE USE OF FFT GRIDS: ================= In case of PAW: --------------- Two FFT grids are used: - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis. It is defined by nfft, ngfft, mgfft, ... Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid. - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres. It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid. In case of norm-conserving: --------------------------- - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ... For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
SOURCE
126 subroutine wfk_analyze(acell, codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, rprim, xred) 127 128 !Arguments ------------------------------------ 129 !scalars 130 character(len=8),intent(in) :: codvsn 131 type(datafiles_type),intent(in) :: dtfil 132 type(dataset_type),intent(in) :: dtset 133 type(pawang_type),intent(inout) :: pawang 134 type(pseudopotential_type),intent(inout) :: psps 135 !arrays 136 real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,dtset%natom) 137 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 138 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 139 140 !Local variables ------------------------------ 141 !scalars 142 integer,parameter :: master=0 143 integer :: comm,nprocs,my_rank,mgfftf,nfftf !,nfftf_tot 144 integer :: optcut,optgr0,optgr1,optgr2,optrad,psp_gencond !,ii 145 !integer :: option,option_test,option_dij,optrhoij 146 integer :: band,ik_ibz,spin,first_band,last_band 147 integer :: ierr,usexcnhat 148 integer :: cplex,cplex_dij,cplex_rhoij,ndij,nspden_rhoij,gnt_option 149 real(dp),parameter :: spinmagntarget=-99.99_dp 150 real(dp) :: ecore,ecut_eff,ecutdg_eff,gsqcutc_eff,gsqcutf_eff,gsqcut_shp 151 !real(dp) :: cpu,wall,gflops 152 !real(dp) :: ex_energy,gsqcutc_eff,gsqcutf_eff,nelect,norm,oldefermi 153 character(len=500) :: msg 154 character(len=fnlen) :: wfk0_path,wfkfull_path 155 logical :: call_pawinit, use_paw_aeur 156 type(hdr_type) :: wfk0_hdr, hdr_kfull 157 type(crystal_t) :: cryst 158 type(ebands_t) :: ebands 159 type(pawfgr_type) :: pawfgr 160 !type(paw_dmft_type) :: paw_dmft 161 type(mpi_type) :: mpi_enreg 162 type(wfd_t) :: wfd 163 type(ddkstore_t) :: ds 164 !arrays 165 integer :: ngfftc(18),ngfftf(18) 166 integer,allocatable :: l_size_atm(:) 167 real(dp),parameter :: k0(3)=zero 168 !real(dp) :: nelect_per_spin(dtset%nsppol),n0(dtset%nsppol) 169 real(dp),pointer :: gs_eigen(:,:,:) 170 complex(gwpc),allocatable :: ur_ae(:) 171 logical,allocatable :: keep_ur(:,:,:),bks_mask(:,:,:) 172 real(dp) :: tsec(2) 173 type(Pawrhoij_type),allocatable :: pawrhoij(:) 174 type(pawfgrtab_type),allocatable :: pawfgrtab(:) 175 !type(paw_ij_type),allocatable :: paw_ij(:) 176 !type(paw_an_type),allocatable :: paw_an(:) 177 type(esymm_t),allocatable :: esymm(:,:) 178 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 179 180 !************************************************************************ 181 182 DBG_ENTER('COLL') 183 184 ! abirules! 185 if (.False.) write(std_out,*)acell,codvsn,rprim,xred 186 187 comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 188 189 wfk0_path = dtfil%fnamewffk 190 if (my_rank == master) then 191 ! Accept WFK file in Fortran or netcdf format. 192 if (nctk_try_fort_or_ncfile(wfk0_path, msg) /= 0) then 193 ABI_ERROR(sjoin("Cannot find GS WFK file:", ch10, msg)) 194 end if 195 end if 196 call xmpi_bcast(wfk0_path, master, comm, ierr) 197 call wrtout(ab_out, sjoin("- Reading GS states from WFK file:", wfk0_path)) 198 199 !call cwtime(cpu,wall,gflops,"start") 200 201 ! Costruct crystal and ebands from the GS WFK file. 202 call wfk_read_eigenvalues(wfk0_path, gs_eigen, wfk0_hdr, comm) !,gs_occ) 203 call wfk0_hdr%vs_dtset(dtset) 204 205 cryst = wfk0_hdr%get_crystal() 206 call cryst%print(header="Crystal structure from WFK file") 207 208 ebands = ebands_from_hdr(wfk0_hdr, maxval(wfk0_hdr%nband), gs_eigen) 209 210 !call ebands_update_occ(ebands, spinmagntarget) 211 call ebands_print(ebands,header="Ground state energies", prtvol=dtset%prtvol) 212 ABI_FREE(gs_eigen) 213 214 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfftc,ngfftf,& 215 gsqcutc_eff=gsqcutc_eff,gsqcutf_eff=gsqcutf_eff,gmet=cryst%gmet,k0=k0) 216 217 call print_ngfft(ngfftc, header='Coarse FFT mesh used for the wavefunctions') 218 call print_ngfft(ngfftf, header='Dense FFT mesh used for densities and potentials') 219 220 ! Fake MPI_type for the sequential part. 221 call initmpi_seq(mpi_enreg) 222 call init_distribfft_seq(mpi_enreg%distribfft,'c',ngfftc(2),ngfftc(3),'all') 223 call init_distribfft_seq(mpi_enreg%distribfft,'f',ngfftf(2),ngfftf(3),'all') 224 225 ! =========================================== 226 ! === Open and read pseudopotential files === 227 ! =========================================== 228 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcutf_eff,pawrad,pawtab,psps,cryst%rprimd,comm_mpi=comm) 229 230 ! ============================ 231 ! ==== PAW initialization ==== 232 ! ============================ 233 if (dtset%usepaw == 1) then 234 call chkpawovlp(cryst%natom,cryst%ntypat,dtset%pawovlp,pawtab,cryst%rmet,cryst%typat,cryst%xred) 235 236 cplex_dij=dtset%nspinor; cplex=1; ndij=1 237 238 ABI_MALLOC(pawrhoij,(cryst%natom)) 239 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,& 240 nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc) 241 call pawrhoij_alloc(pawrhoij,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,cryst%typat,pawtab=pawtab) 242 243 ! Initialize values for several basic arrays 244 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 245 246 ! Test if we have to call pawinit 247 call paw_gencond(dtset,gnt_option,"test",call_pawinit) 248 249 if (psp_gencond==1 .or. call_pawinit) then 250 call timab(553,1,tsec) 251 gsqcut_shp = two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 252 call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,zero,dtset%pawlcutd,dtset%pawlmix,& 253 psps%mpsang,dtset%pawnphi,cryst%nsym,dtset%pawntheta,pawang,Pawrad,& 254 dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero) 255 call timab(553,2,tsec) 256 257 ! Update internal values 258 call paw_gencond(dtset,gnt_option,"save",call_pawinit) 259 260 else 261 if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij =2 262 if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla=2 263 end if 264 265 psps%n1xccc=MAXVAL(pawtab(1:cryst%ntypat)%usetcore) 266 267 ! Initialize optional flags in pawtab to zero 268 ! (Cannot be done in Pawinit since the routine is called only if some pars. are changed) 269 pawtab(:)%has_nabla = 0 270 pawtab(:)%usepawu = 0 271 pawtab(:)%useexexch = 0 272 pawtab(:)%exchmix =zero 273 pawtab(:)%lamb_shielding =zero 274 275 call setsym_ylm(cryst%gprimd,pawang%l_max-1,cryst%nsym,dtset%pawprtvol,cryst%rprimd,cryst%symrec,pawang%zarot) 276 277 ! Initialize and compute data for DFT+U 278 !paw_dmft%use_dmft=dtset%usedmft 279 !call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,& 280 ! .false.,dtset%jpawu,dtset%lexexch,dtset%lpawu,cryst%ntypat,pawang,dtset%pawprtvol,& 281 ! Pawrad,pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu) 282 !ABI_CHECK(paw_dmft%use_dmft==0,"DMFT not available") 283 !call destroy_sc_dmft(paw_dmft) 284 285 if (my_rank == master) call pawtab_print(pawtab, unit=std_out) 286 287 ! Get Pawrhoij from the header of the WFK file. 288 call pawrhoij_copy(wfk0_hdr%pawrhoij,pawrhoij) 289 290 ! Variables/arrays related to the fine FFT grid. 291 ABI_MALLOC(pawfgrtab,(cryst%natom)) 292 call pawtab_get_lsize(pawtab,l_size_atm,cryst%natom,cryst%typat) 293 cplex=1 294 call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat) 295 ABI_FREE(l_size_atm) 296 297 usexcnhat=maxval(pawtab(:)%usexcnhat) 298 ! 0 if Vloc in atomic data is Vbare (Blochl s formulation) 299 ! 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation) 300 call wrtout(std_out,sjoin("using usexcnhat= ",itoa(usexcnhat))) 301 ! 302 ! Identify parts of the rectangular grid where the density has to be calculated === 303 !optcut=0; optgr0=dtset%pawstgylm; optgr1=0; optgr2=0; optrad=1-dtset%pawstgylm 304 !if (dtset%xclevel==2 .and. usexcnhat>0) optgr1=dtset%pawstgylm 305 optcut=1; optgr0=1; optgr1=1; optgr2=1; optrad=1 306 307 call nhatgrid(cryst%atindx1,cryst%gmet,cryst%natom,cryst%natom,cryst%nattyp,ngfftf,cryst%ntypat,& 308 optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,cryst%rprimd,cryst%typat,cryst%ucvol,cryst%xred) 309 310 call pawfgrtab_print(pawfgrtab,cryst%natom,unit=std_out,prtvol=dtset%pawprtvol) 311 312 !ABI_MALLOC(ks_nhat,(nfftf,dtset%nspden)) 313 !ks_nhat=zero 314 else 315 ABI_MALLOC(pawfgrtab,(0)) 316 end if !End of PAW Initialization 317 318 select case (dtset%wfk_task) 319 320 case (WFK_TASK_FULLBZ, WFK_TASK_OPTICS_FULLBZ) 321 ! Read wfk0_path and build WFK in full BZ. 322 if (my_rank == master) then 323 wfkfull_path = dtfil%fnameabo_wfk; if (dtset%iomode == IO_MODE_ETSF) wfkfull_path = nctk_ncify(wfkfull_path) 324 call wfk_tofullbz(wfk0_path, dtset, psps, pawtab, wfkfull_path, hdr_kfull) 325 326 ! Write KB form factors. 327 if (dtset%prtkbff == 1 .and. dtset%iomode == IO_MODE_ETSF .and. dtset%usepaw == 0) then 328 call prtkbff(wfkfull_path, hdr_kfull, psps, dtset%prtvol) 329 end if 330 call hdr_kfull%free() 331 end if 332 call xmpi_barrier(comm) 333 334 if (dtset%wfk_task == WFK_TASK_OPTICS_FULLBZ) then 335 ! Calculate the DDK matrix elements from the WFK file in the full BZ. 336 ! This is needed fo computing non-linear properties in optics as symmetries are not 337 ! implemented correctly. 338 ds%only_diago = .False. 339 call ds%compute_ddk(wfkfull_path, dtfil%filnam_ds(4), dtset, psps, pawtab, ngfftc, comm) 340 call ds%free() 341 end if 342 343 case (WFK_TASK_KPTS_ERANGE) 344 call sigtk_kpts_in_erange(dtset, cryst, ebands, psps, pawtab, dtfil%filnam_ds(4), comm) 345 346 !case (WFK_TASK_LDOS) 347 348 case (WFK_TASK_DDK, WFK_TASK_DDK_DIAGO) 349 ! Calculate the DDK matrix elements from the WFK file 350 ds%only_diago = .False.; if (dtset%wfk_task == WFK_TASK_DDK_DIAGO) ds%only_diago = .True. 351 call ds%compute_ddk(wfk0_path, dtfil%filnam_ds(4), dtset, psps, pawtab, ngfftc, comm) 352 call ds%free() 353 354 case (WFK_TASK_EINTERP) 355 ! Band structure interpolation from eigenvalues computed on the k-mesh. 356 call ebands_interpolate_kpath(ebands, dtset, cryst, [0, 0], dtfil%filnam_ds(4), comm) 357 358 case (WFK_TASK_CHECK_SYMTAB) 359 call wfk_check_symtab(wfk0_path, comm) 360 361 case (WFK_TASK_CLASSIFY) 362 ! Band classification. 363 call read_wfd() 364 365 ABI_MALLOC(esymm,(wfd%nkibz,wfd%nsppol)) 366 use_paw_aeur=.False. ! should pass ngfftf but the dense mesh is not forced to be symmetric 367 368 do spin=1,wfd%nsppol 369 do ik_ibz=1,wfd%nkibz 370 first_band = 1 371 last_band = wfd%nband(ik_ibz,spin) 372 call classify_bands(wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,wfd%ngfft,& 373 cryst,ebands,pawtab,pawrad,pawang,psps,dtset%tolsym,esymm(ik_ibz,spin)) 374 end do 375 end do 376 377 call esymm_free(esymm) 378 ABI_FREE(esymm) 379 380 !case (WFK_TASK_UR) 381 ! ! plot KSS wavefunctions. Change bks_mask to select particular states. 382 ! ABI_MALLOC(bks_mask,(Wfd%mband,Wfd%nkibz,Wfd%nsppol)) 383 ! bks_mask=.False.; bks_mask(1:4,1,1)=.True. 384 ! call wfd%plot_ur(Cryst,Psps,Pawtab,Pawrad,ngfftf,bks_mask) 385 ! ABI_FREE(bks_mask) 386 387 case (WFK_TASK_PAW_AEPSI) 388 ! Compute AE PAW wavefunction in real space on the dense FFT mesh. 389 call read_wfd() 390 391 ABI_CHECK(wfd%usepaw == 1, "Not a PAW run") 392 ABI_MALLOC(paw_onsite, (cryst%natom)) 393 call paw_pwaves_lmn_init(paw_onsite,cryst%natom,cryst%natom,cryst%ntypat,& 394 cryst%rprimd,cryst%xcart,pawtab,pawrad,pawfgrtab) 395 396 ! Use dense FFT mesh 397 call wfd%change_ngfft(cryst,psps,ngfftf) 398 band = 1; spin = 1; ik_ibz = 1 399 400 ABI_MALLOC(ur_ae, (wfd%nfft*wfd%nspinor)) 401 call wfd%paw_get_aeur(band,ik_ibz,spin,cryst,paw_onsite,psps,pawtab,pawfgrtab,ur_ae) 402 ABI_FREE(ur_ae) 403 404 call paw_pwaves_lmn_free(paw_onsite) 405 ABI_FREE(paw_onsite) 406 407 !case ("paw_aeden") 408 409 case default 410 ABI_ERROR(sjoin("Wrong task:", itoa(dtset%wfk_task))) 411 end select 412 413 ! Free memory 414 call cryst%free() 415 call ebands_free(ebands) 416 call wfd%free() 417 call pawfgr_destroy(pawfgr) 418 call destroy_mpi_enreg(mpi_enreg) 419 call wfk0_hdr%free() 420 421 ! Deallocation for PAW. 422 if (dtset%usepaw==1) then 423 call pawrhoij_free(pawrhoij) 424 ABI_FREE(pawrhoij) 425 call pawfgrtab_free(pawfgrtab) 426 !call paw_ij_free(paw_ij) 427 !ABI_FREE(paw_ij) 428 !call paw_an_free(paw_an) 429 !ABI_FREE(paw_an) 430 end if 431 ABI_FREE(pawfgrtab) 432 433 DBG_EXIT('COLL') 434 435 contains
wfk_analyze/read_wfd [ Functions ]
[ Top ] [ wfk_analyze ] [ Functions ]
NAME
read_wfd
FUNCTION
Initialize the wavefunction descriptor
SOURCE
447 subroutine read_wfd() 448 449 ! ************************************************************************* 450 451 ABI_MALLOC(keep_ur, (ebands%mband, ebands%nkpt, ebands%nsppol)) 452 ABI_MALLOC(bks_mask, (ebands%mband, ebands%nkpt, ebands%nsppol)) 453 keep_ur = .False.; bks_mask = .True. 454 455 call wfd_init(wfd,cryst,pawtab,psps,keep_ur,ebands%mband,ebands%nband,ebands%nkpt,dtset%nsppol,bks_mask,& 456 dtset%nspden,dtset%nspinor,ecut_eff,dtset%ecutsm,dtset%dilatmx,wfk0_hdr%istwfk,ebands%kptns,ngfftc,& 457 dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm) 458 459 ABI_FREE(keep_ur) 460 ABI_FREE(bks_mask) 461 462 call wfd%read_wfk(wfk0_path,IO_MODE_MPI) 463 !call wfd%test_ortho(cryst, pawtab) 464 465 end subroutine read_wfd 466 467 end subroutine wfk_analyze